• 数值分析

    基础知识

    实数的浮点表示

    机器数的格式

    浮点数的表示:详见二进制转换成其它格式数据,此处的浮点数精度为64位,即指数部分为11位,尾数部分为52位。标准化的IEEE浮点数表示为:

    ±1.bbb×2p

    对于64位浮点数,精度为εmach=252

    IEEE舍入最近法则

    如果在二进制数右边的第53位是0,则向下舍去(在第52位后面截断)。

    如果第53位是1,则向上进位(在第52位上加1);但是如果在53位的右边的所有已知位都是0,此时当且仅当第52位是1时才在第52位上加1,否则在52位后面截断。

    定义:将IEEE双精度浮点数字记做x,利用舍入最近法则记做fl(x)。

    举例:对于数字9.4,二进制表示为 (1001.0110)2,计算舍入误差。

    将二进制转成浮点数为:1.0010110×23=1.00101100×23,因此第53位是1,又53位的右边也有1,所以需要在52位加1,即加上1×252×23,并且舍弃0.1100×252×23=0.1100×249

    对于0.1100,可以这样计算:

    x=0.110024x=1100.1100(241)x=1100x=0.8

    所以舍入误差为

    fl(x)x=0.8×249+1×252×23=0.2×249

    相对舍入误差

    在IEEE机器算术模型中,fl(x)的相对舍入误差不会比机器精度的一半大:

    |fl(x)x||x|12εmach

    机器数的表示

    举例:fl(9.4)对应的存储在机器中的表达?

    对于指数:3 + 1023 = 1026 = (10000000010)2,再加上符号位0,则为

    (010000000010)2=(402)16

    对于尾数:直接转成16进制可得(2CCCCCCCCCCCD)16

    整体16进制表示为4022CCCCCCCCCCCD。

     

    特殊数的十六进制表示

    机器数例子十六进制数
    +Inf1/07FF0000000000000
    -Inf-1/0FFF0000000000000
    NaN0/0FFFxxxxxxxxxxxxx
    +0 0000000000000000
    -0 8000000000000000

    其中x表示不全为0的位,+0和-0在计算中相同,只有符号位不同。

    当指数全为0时,如0000000000000001表示的数为252×21022=21074,这是最小可表达的数字,在此之下的可表达双精度数字完全不能被表示。

     

    浮点数加法

    对于1和253的加法,对于253,虽然53位是1,但是右边全是0且第52位是0,浮点数表示为0,即

    1+253=1

    举例:计算双精度浮点的求和(1+3×253)-1

    image-20231101144341444

    由于求和后的53位是1,又52位也是1,所以需要在第52位上加上1,即最后结果为:

    image-20231101144503039

    在减去1后,结果为251,因为舍入误差导致最后的结果不为0。

     

    有效数字缺失

    举例:找到二次等式x2+912x=3的两个根。

    使用求根公式:x=b±b24ac2a,负号的值可以找到,但是对于正号的值,存在舍入误差会为0,这就是有效数字缺失(两个近似的数字相减)。可以将正号时的求根公式转成

    x2=b+b24ac2a=(b+b24ac)(bb24ac)2a(bb24ac)=2cbb24ac

    这样便可解得正号的根x2=1.062×1011

     

    求解方程

    二分法

    二分法

    给定初始区间[a,b] 使得f(a)f(b)<0 while (b-a)/2 > TOL c = (a+b) / 2 if f(c) = 0, stop,end if f(a)f(c) < 0 b=c else a=c end end 最终的区间[a,b]中包含一个根 近似根为(a+b)/2

    二分法的求解误差

    |xcr|<ba2n+1

    其中xc=(an+bn)/2是解的最优估计值。

    定义:如果误差小于0.5×10p,解精确到小数点后p位。

    举例:使用二分法求解函数在区间[0,1]上的解,要求解精确到小数点后6位,则至少要迭代多少次。

    ba2n+1=12n+1<0.5×106n20

    至少要20次迭代。

     

    不动点迭代

    函数的不动点

    定义:当g(r)=r,实数r是函数g的不动点。

    不动点迭代

    x0=初始估计

    xi+1=g(xi),其中i=0,1,2,...

    举例:使用不动点迭代求解x3+x1=0的根。

    将方程重写为x=1x3,可以定义g(x)=1x3,还可定义g(x)=1x3=x

    还可以将两边加上2x3,可以得到

    3x3+x=1+2x3x(1+3x2)=1+2x3x=1+2x31+3x2=g(x)

    使用x0=0.5作为初始值,但是第一种没有成功收敛(在0和1之间来回跳跃),第二种和第三种成功收敛并且第三种收敛得更快。

     

    不动点迭代几何

    使用cobweb图描述不动点迭代的几何图示

    image-20231101185232379

    可以看到对于g(x)=1-x^3而言,在迭代时会出现发散的情况,因此最后无法收敛,而是否能够收敛与函数在不动点附近的斜率有关。

     

    不动点迭代的线性收敛

    定义:令ei表示迭代过程中第i步时的误差,如果

    limiei+1ei=S<1

    该方法被称为满足线性收敛,收敛速度为S。

    定理:假设函数g是连续可微函数,g(r)=r,S=|g(r)|<1,则不动点迭代对于一个足够接近r的初始估计,以速度S线性收敛到不动点r。

    证明:令xi表示第i步迭代,根据均值定理,在xi和r之间存在ci,满足

    g(xi)g(r)=g(ci)(xir)

    代入xi+1=g(xi)以及r=g(r)可得

    xi+1r=g(ci)(xir)

    定义ei=|xir|,则有

    ei+1=|g(ci)|ei

    如果S=|g'(r)|小于1,则通过替换g',在r附近有一个小的区间满足|g'(x)|<(S+1)/2,这个值仍然小于1,如果xi恰好出现在该区间,则ci也在该区间,因而

    ei+1S+12ei

    所以,误差以(S+1)/2的速度下降,这意味着limixi=r

     

    定义:如果迭代方法对于一个足够接近r的初值能收敛到r,该迭代方法被称为局部收敛到r。

    举例:解释为什么不动点迭代g(x)=cos(x)收敛

    g'(r)=-sin(r)=-sin(0.74)≈-0.67,绝对值比1小,对应解r≈0.74会把附近的初始猜测值吸引过来。

    举例:使用不动点迭代找出方程 cos x = sin x 的根。

    在方程两边加上x可得

    x+cosxsinx=x

    定义g(x)=x + cos x - sin x,可用不动点迭代法求得x=0.7853981633957986。

     

    精度的极限

    前向和后向误差

    举例:使用二分法计算函数f(x)=x32x2+43x827的根,精确到小数点后6位。

    在使用二分法时,在16步计算后便出现错误,f(0.6666641)=0,无法精确到小数点后6位。即使使用MATLAB提供的fzero.m函数计算仍然无法找到5位以上的正确数位。

     

    定义:假设f是一个函数,r是一个根,意味着满足f(r)=0。假设xa是r的近似值,对于根求解问题,近似xa后向误差|f(xa)|前向误差|rxa|

    方程
    方程求解器

    后向误差在左边或者输入一侧(问题数据),这是我们需要对于问题(函数f)改变的量使得方程平衡并输出近似解xa,这个量是|f(xa)|。前向误差是右边或者输出一侧(问题求解),这是我们对于近似解要做的修正,该误差是|rxa|

     

    在上面的例子中,后向误差接近εmach2.2×1016,前向误差则大约为105,双精度数不能在机器精度的相对误差以下可靠计算,因为后向误差不能被可靠降低,同时前向误差也不能减小。

    方程求解方法的终止条件可以基于前向或者后向误差。还有其他相关的终止条件,诸如计算时间的上限问题的上下文对于我们的选择具有指导作用。

     

    威尔金森多项式

    由于函数在多重根位置上f'为0,所以在重根附近形状很平。正因为如此,分离重根可能会遇到困难,但是多重根问题仅仅是冰山一角,没有多重根问题时也可能也会出现问题,比如威尔金森多项式。

    W(x)=(x1)(x2)(x20)

    根搜索的敏感性

    如果在输人中是一个小误差,f在这种情况下对问题进行求解,造成输出中的大问题,这样的问题被称为敏感性问题。

    假设对输入做了一个小变化εg(x),其中ε很小,令Δr对应根中的变化,则

    f(r+Δr)+εg(r+Δr)=0

    进行泰勒展开得到:

    f(r)+(Δr)f(r)+εg(r)+ε(Δr)g(r)+O((Δr)2)=0

    忽略高阶项,得到

    (Δr)(f(r)+εg(r))f(r)εg(r)=εg(r)

    或者,由于ε很小,所以εg'(r)很小。

    Δrεg(r)f(r)+εg(r)εg(r)f(r)

    根的敏感公式

    假设r是函数f(x)的根,并且r+Δr是f(x)+εg(r)的根,则当ε<<f'(r)时,

    Δrεg(r)f(r)

    举例:估计P(x)=(x1)(x2)(x3)(x4)(x5)(x6)106x7的最大根。

    f(x)=(x1)(x2)(x3)(x4)(x5)(x6)ε=106g(x)=x7。如果没有εg(x),则有最大根r=6,但当加上这一项后,根变成了

    r+Δr=6ε675!=62332.8ε=6.0023328

    误差放大因子 = 相对前向误差 / 相对后向误差

     

    条件数也是误差放大度量的一种方式,数值分析是对算法的研究,算法把定义问题的数据作为输入,对应的结果作为输出。条件数指的是理论问题本身所带来的误差放大部分,和用于求解问题的特定算法无关。注意到条件数仅仅度量由于问题本身带来的误差放大,这点很重要。和条件一起,还有一个平行概念,即稳定。稳定指的是由于算法小的输入误差造成的放大,而不是问题本身。如果一个算法在小的后向误差存在的时候,总能给出一个近似解,则称该算法是稳定的。如果问题的条件好,算法稳定,我们可以期望同时具有小的后向误差和前向误差。

    前面的误差放大例子表明根求解过程对于特定的输入变化的敏感性。问题可能或多或少地敏感,依赖于如何设计输入的变化。问题的条件数定义为所有输入变化,或者至少规定类型的变化所造成的最大误差放大。条件数高的问题称为病态问题,条件数在1附近的问题称为良态问题。

     

    牛顿方法

    牛顿方法

    x0=初始估计

    xi+1=xif(xi)f(xi),i=0,1,2,...

     

    定义:令ei表示一个迭代方法第i步后得到的误差,如果满足下式,则该迭代为二次迭代

    M=limiei+1ei2<

    定理:令f是二阶连续可微函数,f(r)=0。如果f'(r)≠0,则牛顿方法局部二次收敛到r,第i步的误差ei满足

    limiei+1ei2=M

    证明参见原书,注意此处

    M=|f(r)2f(r)|

    注意,该定理并不能保证牛顿方法每次都能二次收敛。如求f(x)=x2的实根,此时牛顿迭代公式为

    xi+1=xi2

    误差公式是ei+1=ei/2,是线性收敛。

    定理:假设在区间[a,b]上,(m+1)阶连续可微函数f在r点有一个m阶多重根,则牛顿方法局部收敛到r,第i步误差ei满足

    limiei+1ei=S=m1m

    定理:如果在[a,b]区间上f是(m+1)阶连续函数,包含m>1的多重根,则改进的牛顿方法为:

    xi+1=ximf(xi)f(xi)

    收敛到r,并具有二次收敛速度。

    牛顿方法如同FPI(不动点迭代),可能不会收敛到根.下面的例子是一种可能不收敛的情况。

    举例:对函数f(x)=4x46x211/4使用牛顿方法,初始估计x0=1/2

     

    不需要导数的根求解

    除了重根,牛顿方法比二分法和FPI方法的收敛速度更快。它达到了这种更快的速度是因为使用了更多的信息,尤其是通过函数导数得到的函数切线方向的信息。在某些情况下,可能难以计算导数。在这种情况下,割线方法是牛顿方法的一个非常好的替代。它使用近似值割线替代了切线,并且收敛速度差不多快。

    割线方法及其变体

    xi处的导数的近似可写为差商:

    f(xi)f(xi1)xixi1

    割线方法

    x0x1=初始估计

    xi+1=xif(xi)(xixi1)f(xi)f(xi1),i=1,2,3,

    割线方法以超线性的速度收敛到一个单根,意味着它在线性和二次收敛方法之间。

     

    试位方法和二分法相似,将中点从a和b的平均值改为割线方法定义的点,注意f(a)f(b)<0

    c=af(a)(ab)f(a)f(b)=bf(a)af(b)f(a)f(b)

    c对应的便是a和b所定义的直线与x轴的交点。

    试位方法开始表现得比二分法和割线方法都要好,具有二者最好的性质。但是,二分法在每一步中可以确保消除1/2的不确定性,试位方法却没有能力做出这样的保证,而且在一些例子中可能收敛很慢。

     

    Muller方法是割线方法在不同方向的推广。该方法不是计算经过先前两个点的直线和x轴的交点,而是使用三个前面生成的点x0x1x2,画出通过它们的抛物线y=p(x),并计算抛物线和x轴的交点。一般来讲,抛物线会生成0个或者2个交点。如果有两个交点,接近上一步中的x2的点会被选作x3,通过简单的二次公式计算,就可以确定两种可能。如果抛物线和x轴不相交,就会出现复数解,能够处理复数代数的软件就可以计算对应的解。

    逆二次插值(IQI)是割线方法到抛物线的一种相近的泛化方法,但是使用形如x=p(y)的抛物线,而不是 Muller方法所使用的y=p(x),这个抛物线和x轴只有一个交点。经过三点(a,A),(b,B),(c,C)的二阶多项式x=P(y)为:

    P(y)=a(yB)(yC)(AB)(AC)+b(yA)(yC)(BA)(BC)+c(yA)(yB)(CA)(CB)

    其中P(A)=a,P(B)=b,P(C)=c,用y=0带入可得得到抛物线和x轴的交点,可以得到

    P(0)=cr(rq)(cb)+(1r)s(ca)(q1)(r1)(s1)

    其中q=f(a)/f(b),r=f(c)/f(b),s=f(c)/f(a)。

    逆二次插值

    a=xib=xi+1c=xi+2A=f(xi)B=f(xi+1)C=f(xi+2)

    下一步的xi+3=P(0)

    xi+3=xi+2r(rq)(xi+2xi+1)+(1r)s(xi+2xi)(q1)(r1)(s1)

    可以用新的xi+3替代最旧的估计xi,另一种方式使用最新估计替换最近的三个估计中的后向误差最大的那个。

    Brent方法

    用于连续函数f,区间的边界是a和b,同时f(a)f(b)<0。Brent方法记录当前点xi,该点具有最优的后向误差,同时有包含根的区间[ai,bi]。简单来讲,尝试使用逆二次方法,并在下述情况下,使用结果来替代xiaibi中的一个:

    (1)后向误差得到改进;

    (2)包含根的区间至少减小一半;

    否则,尝试使用割线方法以实现相同的目的。如果割线方法也失败了,则使用二分法,保证至少减少一半的不确定性。

     

    方程组

    高斯消去法不过多介绍

    LU分解

    可以将高斯消去法的消去步骤表示为LU分解,其中L是单位矩阵进行高斯消去操作后的下三角矩阵,U是原矩阵经过高斯消去得到的矩阵。

    一旦知道L和U,问题Ax=b可以写作LUx=b,定义新的“辅助”向量c=Ux,则回代是个两步的过程

    (1)对于方程Lc=b,求解c

    (2)对于方程Ux=c,求解x

    但并不是所有的矩阵都能进行LU分解,如

    [0111]

    无法被LU分解。

     

    误差来源

    在高斯消去法中有两个潜在的误差来源,第一个是病态问题的概念和解对于输入数据的敏感性相关,第二个是淹没,在大部分问题中可以通过一个简单的修正,称为部分主元,进行避免。

    误差放大和条件数

    定义:向量x=(x1,,xn)的无穷范数或者最大范数为x=max|xi|,i=1,...,n,即x所有元素的最大值。

    定义:令xa是线性方程组Ax=b的近似解,余项是向量r=b-Axa后向误差是余项的范数bAxa前向误差xxa

    举例:找出近似解[-1,3.0001]的后向误差和前向误差,方程组如下:

    x1+x2=21.0001x1+x2=2.0001

    精确解为[1,1],后向误差计算:

    bAxa=[22.0001][111.00011][13.0001]=[0.00010.0001]

    后向误差为0.0001,前向误差计算:

    xxa=[11][13.0001]=[22.0001]

    前向误差为2.0001。可以看到前向误差和后向误差差别很大。下面这张图可以帮助理解这一情况

    image-20231102104200536

    可以看到(1,1)是它们的交点,不过(-1,3.0001)也差一点在两条直线上(图中的差异是人为放大了1000倍,实际上离得非常近)。

    把余项表示为r=bAxa,系统Ax=b的相对后向误差定义为

    rb

    相对前向误差定义为

    xxax

    误差放大因子则是二者的比值,即 相对前向误差/相对后向误差。

    定义:方阵A的条件数cond(A)为求解Ax=b时,对于所有右侧向量b,可能出现的最大误差放大因子。

    方阵A的矩阵范数定义为每行元素绝对值求和,其中的最大值作为矩阵A的范数。

    定理:n×n矩阵A的条件数是:

    cond(A)=AA1

    举例:计算矩阵A的条件数,其中

    A=[111.00011]

    计算A=2.0001,A的逆为

    A1=[10000100001000110000]

    A1=20001,A的条件数是cond(A)=40004.0001,在这个系统中,对于任何其他的b,误差放大因子将小于或者等于40004.0001。

    在浮点算术中,相对后向误差不可能小于εmach,这是由于b的元素的存储已经引入了和εmach差不多大的误差,在求解Ax=b可能出现的相关前向误差是εmach·cond(A)。换句话讲,如果cond(A)10k,我们在计算x时,将丢掉k位数字精度。

    以上例为例,cond(A)4×104,对于双精度16位精确度而言,只能求解到大约12位的精确数字(实际上会更低,matlab只能求解到11位精确值)。

    对于希尔伯特矩阵H(Hij=1/(i+j1)),其对应的条件数非常大。

    对于向量和矩阵范数,有A+BA+B

    对于向量x的1-范数为x1=|x1|++|xn|,n×n矩阵A的矩阵1-范数是A1的最大绝对列和,即列向量的1-范数的最大值。

     

    淹没

    举例:考虑方程组

    1020x1+x2=1x1+2x2=4

    (1)精确解约为[2,1]

    (2)IEEE双精度,

    [10201|112|4]r21020r1[10201|1021020|41020]

    由于存在舍入,2102041020存为1020,因此解为[0,1]

    (3)IEEE双精度,使用行交换

    [12|410201|1]r21020r1[12|4012×1020|14×1020]

    同样存在舍入导致12×102014×1020存为1,解得[2,1],注意这里的结果是确切的结果,不像精确解存在微小的差距。

    第2种方法在消去过程中的乘子1020,从底部的方程减去顶部方程的1020倍,底部方程会被抑制,或者称为"swamp,",尽管初始的时候有两个独立的方程或者源信息,但是经过第2种计算的消去,这里仅仅留下 了顶部方程的两个副本。由于底部方程消失,出于所有实际的目的,我们不能指望计算结果可以满足底部的方程,实际上得到的解也不能满足对应的方程。

    而另一方面,第3种方法在消去过程中没有产生覆盖,因为乘子是1020, 在消去之后,原始的两个方程仍然存在,并且变为三角形式,结果是更加精确的近似解。

    高斯消去的过程中保证乘子尽可能小,同时避免淹没。幸运的是,通过对朴素的高斯消去法的简单修正,可以使得高斯消去法中的乘子的绝对值不大于1。这种新的原则,包括在表格中明智的行交换,称为部分主元方法,在下一节中将进行详细的描述。

     

    PA=LU分解

    部分主元

    部分主元包含在每一步消去步骤之前的比较,找到第一列中最大的一个元素,其对应行和主元行进行交换,在当前情况下,主元行是第一行。在消去过程中,对每一列都实施相同的策略。在消去第k列时,找到第p行,k≤p≤n,定位最大的|apk|,必要时交换第p行和第k行,然后继续进行消去。注意到使用部分主元方法保证所有乘子,或者L的元素的绝对值不大于1。通过这种对于高斯消去方法所进行的小的改变,淹没问题可以完全避免。

    举例:使用高斯消去的部分主元方法求解方程组

    x1x2+3x3=3x12x3=12x1+2x2+4x3=0

    由于第一列中a31最大,所以交换第一行和第三行

    [113|3102|1224|0]r1r3[224|0102|1113|3]r2+12r1[224|0010|1113|3]r312r1[224|0010|1021|3]

    在消去第二列时,由于a32a22大,所以交换第二行和第三行

    [224|0010|1021|3]r2r3[224|0021|3010|1]r3+12r2[224|0021|3000.5|0.5]

    解得x=[1,1,-1]

     

    置换矩阵

    定义:置换矩阵是一个n×n的矩阵,其在每一行、每一列仅有一个1,其他全部为0。

    定理:令P是通过对单位矩阵实施一组特定的行交换后得到的一个n×n的置换矩阵。则对于任意的n×n矩阵 A,PA对应于对矩阵A实施同样的行交换得到的结果。

     

    PA=LU分解

    我们把关于高斯消去所知道的一切组合起来,进行PA=LU分解。这是部分主元的高斯消去的矩阵形式,PA=LU分解是求解线性方程组的主要方法。正如名字中所暗示的,PA=LU是对于矩阵A包含行交换的LU分解。在部分主元中,开始的时候我们并不知道需要进行置换的列,所以我们必须谨慎地把行置换的信息加 入分解中,我们需要记录高斯消去法中所有的主元。

    使用PA=LU求解方程组Ax=b时

    PAx=PbLUx=PbLc=PbcUx=cx

    其中P即为单位阵进行部分主元法交换行后产生的置换矩阵。

     

    迭代方法

    雅可比方法

    雅可比(Jacobi)方法是方程组系统中的一种形式的不动点迭代。在FPI中,第一步是重写方程,进而求解未知量。雅可比方法以如下标准方式进行该重写步骤:求解第i个方程得到第i个未知变量,然后使用不动点迭代,从初始估计开始,进行迭代。

    举例:使用雅可比方法求解下面的线性方程组

    3x1x2+2x3=5x1+4x22x3=42x1+x2+4x3=0

    求解未知变量得到

    x1=x22x353x2=x1+2x3+44x3=2x1x24

    迭代过程为:

    [x1x2x3]=[000][x1x2x3]=[x22x353x1+2x3+442x1x24]=[534404]=[5310][x1x2x3]=[x22x353x1+2x3+442x1x24]=[12×05353+2×0+442×5314]=[43712712]...[x1x2x3]=[1.8330.8890.694]

    定义:当对角线上的元素的绝对值比该行上其它元素的绝对值之和还要大的时候,认为n×n的矩阵A=(aij)严格对角占优矩阵。

    定理:如果n×n矩阵A是严格对角占优矩阵,则(1) A是非奇异矩阵,(2)对于所有向量b和初始估计,对Ax=b应用雅可比方法都会收敛到(唯一)解。

    但是这个只是充分条件,不代表不满足严格对角占优矩阵,就不会收敛到唯一解,上例中就不是严格对角占优矩阵,可以看到第一行的对角线元素为3,只是等于第一行其它元素的绝对值之和。

    雅可比方法是不动点迭代的一种形式,令D表示A的主对角线矩阵,L表示矩阵A的下三角矩阵(主对角线以下的元素),U表示上三角矩阵(主对角线以上的元素),则A=L+D+U,求解的方程变为Lx+Dx+Ux=b。 注意,这里对于L和U的使用和LU分解中不同,当前L和U的主对角线元素都是零。

    Ax=b

    (D+L+U)x=b

    Dx=b-(L+U)x

    x=D1(b-(L+U)x)

    令g(x)=D1(b-(L+U)x)

    雅可比方法

    x0=初始向量

    xk+1=D1(b-(L+U)xk), k = 0,1,2,...

    高斯-塞德尔方法和SOR

    高斯-塞德尔方法和雅可比方法的唯一差异是在前者中,最近更新的未知变量的值在每一步中都使用,即 使更新发生在当前步骤。在使用相同步骤的情况下,该近似更加精确。如果收敛,高斯-塞德尔方法常常比雅可比方法收敛更快。

    高斯-塞德尔方法

    x0=初始向量

    xk+1=D1(b-Uxk-Lxk+1) k = 0,1,2,...

    举例:对如下系统应用高斯-塞德尔方法

    3uv+2w=5u+4v2w=42u+v+4w=0

    迭代公式为

    [uk+1vk+1wk+1]=[vk2wk53uk+1+2wk+442uk+1vk+14]

    连续过松弛(SOR)方法使用高斯-塞德尔的求解方向,并使用过松弛以加快收敛速度。令w是一个实数,并将新的估计中的每个元素xk+1定义为w乘上高斯-塞德尔公式和1-w乘上当前估计xk的平均。数字w被称为松弛参数,而当w>1时被称为过松弛。

    使用SOR后的迭代公式

    [uk+1vk+1wk+1]=[(1w)uk+wvk2wk53(1w)vk+wuk+1+2wk+44(1w)wk+w2uk+1vk+14]

    可以将系统看成不动点迭代,如下:

    Ax=b(wL+wD+wU)x=wb(wL+D)x=wbwUx+(1w)Dxx=(wL+D)1[(1w)DxwUx]+w(wL+D)1b

    连续过松弛(SOR)

    x0=初始向量

    xk+1=(wL+D)1[(1w)DxkwUxk]+w(wL+D)1b

     

    收敛性的证明略

     

    稀疏矩阵计算

    基于高斯消去的直接方法,通过有限步数的计算就可以得到解。为什么我们要研究并使用迭代方法,特别是迭代方法只能得到近似解,并需要多步计算以得到收敛。有两个主要的原因让我们使用诸如高斯-塞德尔这样的迭代方法。两个原因都是基于如下的事实:迭代方法中的一步计算,仅仅需要LU分解的浮点计算所需要时间的一小部分。对于n×n的矩阵做一次高斯消去需要n3次的操作,雅可比方法的一步,仅需要大约n2次的乘法(对于每个矩阵元素做一次乘法操作)以及大约相同数量的加法。

    如果已知解的较好的近似,在这种特定的情况下支持使用迭代技术。例如,假设知道Ax=b的解,随后A与b同时或者单个仅仅发生小的变化。我们可以想象一个动态系统,当A和b改变时,对A和b进行待续的度量,因而需要待续更新精确的解向量x。如果把前面问题的解作为新的相似问题的初始估计,使用雅可 比或者高斯-塞德尔可以得到更快的收敛。这种技术通常称为修饰,这是由于该问题从一个近似解开始,该近似解可能对应前面一个相关问题的解,然后仅仅修饰近似解使其更加精确。

    第二个使用迭代方法的主要原因是用于求解稀疏的方程组。当矩阵中的很多元素都是0系数矩阵被称为稀 疏矩阵。通常,对于稀疏矩阵中的n2个矩阵元素,只有O(n)个非零元素。一个完全的矩阵恰恰相反,其中没有几个元素可以假设为0。对稀疏矩阵使用高斯消去经常会导致填充,这使得其中系数矩阵由于必要的行交换从稀疏变为不稀疏。由于这个原因,高斯消去以及PA=LU实现的效率对于稀疏矩阵变得可疑,在这种情况下迭代方法是一种合理的选择。

     

    用来对称正定矩阵的方法

    对称正定矩阵

    定义:如果AT=A,则n×n矩阵A是对称矩阵,如果对于所有向量x≠0,xTAx>0,则矩阵A是正定矩阵

    性质:如果n×n矩阵A是对称矩阵,则A是正定矩阵,当且仅当所有特征值是正数。

    性质:如果A是n×n对称正定矩阵,X是一个满秩n×m矩阵,n≥m,则XTAX是m×m对称正定矩阵。

    定义:方阵A的主子矩阵是一个方的子矩阵,其对角线元素是矩阵A的对角线元素。

    性质:任何对称正定矩阵的主子矩阵也是对称正定矩阵。

     

    Cholesky分解

    考虑一个对称正定矩阵

    [abbc]

    将A写成RTR的形式

    [abbc]=[a0uv][au0v]=[auauau2+v2]

    与上式作比较,u=b/a以及v2=cu2,Cholesky分解如下:

    A=[abbc]=[a0bacb2/a][aba0cb2/a]=RTR

    定理:如果A是n×n对称正定矩阵,则存在上三角n×n矩阵R满足A=RTR

    Cholesky分解

    for k = 1,2,...,n

    if Akk<0,stop,end

    Rkk=Akk

    uT=1RkkAk,k+1:n

    Rk,k+1:n = uT

    Ak+1:n,k+1:n=Ak+1:n,k+1:nuuT

    end

    举例:计算如下矩阵的Cholesky分解

    [4222242411]

    k=1:R11=A11=2uT=A1,2:3/R11=[2,2]/2=[1,1]

    R1,2:3=[1,1],对A2:3,2:3进行运算:

    [24411][1111]=[13310]

    k=2:现在在这个2×2的子矩阵中重复这个过程,找到R22=A22=1以及uT=A23=3R23=uT=3,下面的1×1主子矩阵是10-(-3)(-3)=1,因而R33=1,所以R等于

    R=[211013001]

    共轭梯度方法

    共扼梯度的思路依赖于内积思想的推广,因为(v,w)=(w,v),以及对于标量α和β,有(αv+βw,u)=α(v,u)+β(w,u),所以欧几里得内积(v,w)=vTw对称并对于输入v和w线性,欧几里得内积也是正定的,当v≠0时,(v,v)>0。

    定义:令A是对称正定的n×n矩阵,对于两个n维向量v和w,定义A内积

    (v,w)A=vTAw

    (v,w)A=0时,向量v和w为A共轭

    共轭梯度方法

    x0=初始估计

    d0=r0=bAx0

    for k = 0,1,2,...,n-1

    if rk=0,stop,end

    αk=rkTrkdkTAdk

    xk+1=xk+αkdk

    rk+1=rkαkAdk

    βk=rk+1Trk+1rkTrk

    dk+1=rk+1+βkdk

    end

    xk表示当前解,rk表示余项,dk表示更新xk得到改进的xk+1时所使用的新的搜索方向。该方法能够成功的原因在于所有的余项都和前面的余项正交。如果能做到所有的余项正交,该方法搜索所有的正交方向,在经过至多n步就可以得到余项为零的正确解。实现所有余项的正交的关键在于选择搜索方向dk,并使之两两共扼。共轭的概念推广了正交的概念,并据此在算法的名字中也包含“共轭”。

    对于αkβk的选择,是为了保证下一余项向量和前面所有的余项向量都正交。选择αk使得新的余项rk+1和方向dk正交:

    dkTrk+1=dkTrkαkdkTAdk=0αk=rkTrkdkTAdk

    对于βk使dk1rk正交,有

    dkrk=βk1dk1rkTdkrkTrk=0

    rkTdk=rkTrk,为了保证dk两两A共轭,

    dk+1=rk+1+βkdkdkTAdk+1=dkTArk+1+βkdkTAdk0=dkTArk+1+βkdkTAdkβk=dkTArk+1dkTAdk

    βk还可以继续化简成伪代码中的形式,具体的证明参考原书。

     

    预条件

    通过使用预条件技术,可以使得诸如共扼梯度方法的迭代方法收敛速度加快。迭代方法的收敛率通常直接或者间接依赖于系数矩阵A的条件数。预条件方法的思想是降低问题中的条件数。

    n×n的线性方程组Ax=b的预条件形式是

    M1Ax=M1b

    其中M是可逆的n×n矩阵,称为预条件子。我们所要做的就是在方程两侧左乘上该矩阵,预条件矩阵试图对矩阵A逆转,从而可以有效地降低问题的条件数。概念上来讲,它试图同时做两件事:矩阵M应该(1)和矩阵A足够接近,(2)容易求逆。这两个目的常常彼此对立。

    一种特别简单的方式是雅可比预条件子M=D,其中D是A的对角线矩阵。D的逆矩阵是D的元素的倒数.例如在一个严格对角占优矩阵中,雅可比预条件子和A相似,同时非常容易求逆。注意到,对称正定矩阵的每个对角线元素都严格为正,因而计算倒数不是问题。

    预条件共轭梯度法

    x0=初始估计

    r0=bA0

    d0=z0=M1r0

    for k = 0,1,2,...,n-1

    if rk=0,stop,end

    ak=rkTzk/dkTAdk

    xk+1=xk+αkdk

    rk+1=rkαkAdk

    zk+1=M1rk+1

    βk=rk+1Tzk+1/rkTzk

    dk+1=zk+1+βkdk

    end

    在对称连续过松弛(SSOR)中,预条件子定义如下:

    M=(D+wL)D1(D+wU)

    其中A=L+D+U被分为下三角部分、对角线以及上三角部分。w是一个在0和2之间的常数。在特例中w=1,这被称为高斯-塞德尔预条件子

    注意到SSOR预条件子定义为下三角矩阵和上三角矩阵的乘积M=(I+wLD1)(D+wU),因而方程z=M1v可以通过两次回代求解。

    (I+wLD1)c=v(D+wU)z=c

    非线性方程组

    多元牛顿方法

    多变量情况下函数导数f'对应的是雅可比矩阵,定义为

    DF(x)=[f1uf1vf1wf2uf2vf2wf3uf3vf3w]

    多变量牛顿方法

    x0=初始向量

    xk+1=xk(DF(xk))1F(xk),k=0,1,2,...

    由于出现计算矩阵的逆,所以令xk+1=xks,其中s是DF(xk)s=F(xk)的解,所以迭代步骤为

    {DF(xk)s=F(xk)xk+1=xk+s

    Broyden方法

    如果可以计算雅可比矩阵,那么牛顿方法是一个好选择,那么最好的替代方法就是Broyden方法。

    Broyden方法

    x0=初始向量

    A0=初始矩阵

    for i = 0,1,2,...

    xi+1=xiAiF(xi)

    Ai+1=Ai+(Δi+1Aiδi+1)δi+1Tδi+1Tδi+1

    end

    其中δi+1=xi+1xi,Δi+1=F(xi+1)F(xi)

    Broyden方法Ⅱ

    x0=初始向量

    A0=初始矩阵

    for i = 0,1,2,...

    xi+1=xiBiF(xi)

    Bi+1=Bi+(δi+1BiΔi+1)δi+1TBiδi+1TBiδi+1

    end

    其中δi+1=xi+1xi,Δi+1=F(xi+1)F(xi)

    首先,需要初始向量x0和初始估计B0,如果难以计算导数,可以使用B0=I。Broyden方法 Ⅱ的一个可以察觉的缺点是在一些应用中需要估计雅可比矩阵,但是这个矩阵并不容易得到。矩阵Bi是对雅可比矩阵逆的估计。而Broyden方法I正相反,一直记录了用来Ai估计雅可比。两个版本的Broyden方法都超线性收敛到单根,比牛顿方法的二次收敛要慢一些。

     

    插值

    数据和插值函数

    定义:如果对于每个1≤i≤n,P(xi)=yi,则称函数y=P(x)插值数据点(x1,y1),...,(xn,yn)

    拉格朗日插值

    假设给出三点(x1,y1),(x2,y2),(x3,y3),则多项式

    P2(x)=y1(xx2)(xx3)(x1x2)(x1x3)+y2(xx1)(xx3)(x2x1)(x2x3)+y3(xx1)(xx2)(x3x1)(x3x2)

    一般来说,如果有n个点(x1,y1),...,(xn,yn),对于1~n之间的每个k,定义n-1次多项式

    Lk(x)=(xx1)(xxk1)(xxk+1)(xxn)(xkx1)(xkxk1)(xkxk+1)(xkxn)

    Lk(x)具有一个有趣的性质,即Lk(xk)=1,而Lk(xj)=0jk。定义了n-1次多项式

    Pn1(x)=y1L1(x)+ynLn(x)

    定理(多项式插值的主定理):令(x1,y1),...,(xn,yn)是平面中的n个点,具有不同的xi坐标,则存在一个并且仅有一个n-1次或者更低次的多项式P满足P(xi)=yi,i=1,...,n。

     

    牛顿差商

    定义:用f[x1...xn] 表示(唯一)多项式的xn1项的系数,该多项式插值(x1,f(x1)),...,(xn,f(xn))。

    牛顿差商公式

    P(x)=f[x1]+f[x1x2](xx1)+f[x1x2x3](xx1)(xx2)+f[x1x2x3x4](xx1)(xx2)(xx3)+...+f[x1x4](xx1)(xxn1)

    牛顿差商

    给定x=[x1,...,xn],y=[y1,...,yn]

    for j = 1,...,n

    f[xj]=yj

    end

    for i = 2,...,n

    for j=1,...,n+1-i

    f[xj...xj+i1]=(f[xj+1...xj+i1]-f[xj...xj+i2])/(xj+i1-xj)

    end

    end

     

    插值代码

     

    插值误差

    不过多赘述。

     

    切比雪夫插值

    切比雪夫插值是一种特定最优的点间距选取方式。

    切比雪夫理论

    切比雪夫插值的动机是在插值区间上,提高对如下插值误差的最大值的控制

    (xx1)(xx2)(xxn)n!f(n)(c)

    从现在开始让我们把区间固定在[-1,1]。多项式插值误差的分子

    (xx1)(xx2)(xxn)

    本身是一个关于x的n阶多项式,并在区间[-1,1]上具有极值。是否可能在区间[-1,1]找到特定的x1,...,xn使得多项式插值误差的分子的最大值足够小?这被称为插值误差的最小最大问题。

    定理:选择实数1x1,,xn1,使得

    max1x1|(xx1)(xxn)|

    尽可能小,则

    xi=cos(2i1)π2n,i=1,n

    对应的最小值是1/2n1,实际上,通过

    (xx1)(xxn)=12n1Tn(x)

    可以得到极小值,其中Tn(x)表示n阶切比雪夫多项式。

    切比雪夫多项式

    定义n阶切比雪夫多项式Tn(x)=cos(narccosx),对于每个n都是关于x的多项式。令y=arccosx,因而cosy=x。

    Tn+1(x)=cos(n+1)y=cos(ny+y)=cosnycosysinnysinyTn1(x)=cos(n1)y=cos(nyy)=cosnycosy+sinnysinyTn+1(x)+Tn1(x)=2cosnycosy=2xTn(x)Tn+1(x)=2xTn(x)Tn1(x)

    事实deg(Tn)=n,主导系数是2n1

    事实Tn(1)=1Tn(1)=(1)n

    事实Tn(x)的最大绝对值是1,-1≤x≤1,这由Tn(x)=cosy的形式得到。

    事实Tn(x)的所有过零点都在-1~1之间,实际上,过零点是0=cos(n arccosx)的解,由于cosy=0,当且仅当y=奇数∙(π/2)。

    narccosx=oddπ/2x=cosoddπ2n

    事实Tn(x)在-1和1之间一共往返变化n+1次,实际上,这发生在cos0,cosπ/n,...,cos(n-1)π/n,cosπ。

     

    区间的变化

    关于切比雪夫插值的讨论局限在区间[-1,1],可以将方法推广到一般的区间[a,b],移动基点使得它们在区间[a,b]上的相对位置和在区间[-1,1]上一致。这可以通过如下两步实现:

    (1)使用因子(b-a)/2拉伸点(这是两个区间长度的比值),

    (2)将点平移(b+a)/2,使得中心从0移动到区间[a,b]的中心。

    换句话讲,从原始点

    cosoddπ2nba2cosoddπ2n+ba2

    使用区间[a,b]上新的切比雪夫基点x1,...,xn,插值误差公式中分子部分的上界也发生了改变,这是由于对因子xxi拉伸(b-a)/2,结果最小最大值1/2n1必须被替换为[(ba)/2]n/2n1

    切比雪夫插值节点

    在区间[a,b]

    xi=b+a2+ba2cos(2i1)π2n

    i=1,...,n,不等式

    |(xx1)(xxn)|(ba2)n2n1

    在区间[a,b]上成立。

    举例:在区间[0,π/2]找出4个切比雪夫基点进行插值,找出切比雪夫插值误差的上界,在区间中f(x)=sinx。

    xi=π2+02+π202cos(2i1)π8,i=1,2,3,4
    x1=π4+π4cosπ8x2=π4+π4cos3π8x3=π4+π4cos5π8x4=π4+π4cos7π8

    插值误差的最坏情况是

    |sinxP3(x)|=|(xx1)(xx2)(xx3)(xx4)|4!|f(c)|(π202)44!2310.00198

    这里获得的只是插值时使用的基点,即牛顿差商x1x2,...。

     

    三次样条

    样条是另一种数据插值的方式,在多项式插值中,多项式给出的单一公式满足所有数据点,而样条 则使用多个公式,其中每个都是低阶多项式,来通过所有数据点。线性样条可以成功地对任意的n 个点集进行插值,但是线性样条函数缺乏平滑,三次样条则可以解决线性样条的这个缺点。三次样条在两个数据点之间使用3阶(cubic)多项式替换线性样条两点之间的线性函数。由于通过n个点的三次样条无穷多,所以需要添加额外条件。

    样条的性质

    假设给定n个点(x1,y1),...,(xn,yn),其中xi不同且升序。通过n个点的三次样条S(x)是一组三次多项式

    S1(x)=y1+b1(xx1)+c1(xx1)2+d1(xx1)3x[x1,x2]S2(x)=y2+b2(xx2)+c2(xx2)2+d2(xx2)3x[x2,x3]Sn1(x)=yn1+bn1(xxn1)+cn1(xxn1)2+dn1(xxn1)3x[xn1,xn]

    并具有如下性质:

    (1)Si(xi)=yi,Si(xi+1)=yi+1,i=1,2,...,n-1。

    (2)Si1(xi)=Si(xi),其中i=2,...,n-1。

    (3)Si1(xi)=Si(xi),其中i=2,...,n-1。

    性质1保证样条S(x)插值数据点;性质2使得相邻的样条段在它们相遇的地方斜率相同;性质3则保证在两条样条段相邻的地方曲率相同,该曲率由二阶导数表示。

    只有满足上面的这些条件才被称为三次样条

    性质4a(自然样条)S1(x1)=Sn1(xn)=0

    满足性质4a这两个附加条件的三次样条被称为自然三次样条,自然三次样条是唯一的。

    δi=xi+1xiΔi=yi+1yi

    image-20231103222300906

    自然三次样条

    给定x=[x1,,xn],其中x1<<xny=[y1,,yn]

    for i = 1,...,n-1

    ai=yi

    δi=xi+1xi

    Δi=yi+1yi

    end

    求解3.24得到c1,,cn

    for i = 1,...,n-1

    di=ci+1ci3δi

    bi=Δiδiδi3(2ci+ci+1)

    end

    除了性质4a给出的自然三次样条,还有曲率调整三次样条、钳制三次样条、抛物线端点的三次样条、非纽结三次样条。

     

    贝塞尔曲线

    贝塞尔样条是一个允许用户控制节点处斜率的样条。作为额外自由控制的代价,在节点处的一阶和二阶导数的平滑性不再能保证,而这种平滑性是前面三次样条本身就具有的性质。贝塞尔样条适合不时出现角点(一阶导数不连续)和曲率突变(二阶导数不连续)的情况。

    平面贝塞尔样条的每一段由4个点(x1,y1)(x2,y2)(x3,y3)(x4,y4)所确定,第一个点和最后一个点是样条的起点和终点,中间的两个点是控制点。如图所示,曲线以(x2x1,y2y1)离开(x1,y1),并以切线方向(x4x3,y4y3)(x4,y4)点结束。

    image-20231103223127691

    贝塞尔曲线

    给定端点(x1,y1)(x4,y4)

    控制点(x2,y2)(x3,y3)

    bx=3(x2x1)cx=3(x3x2)bxdx=x4x1bxcxby=3(y2y1)cy=3(y3y2)bydy=y4y1bycy

    定义在0≤t≤1的贝塞尔曲线如下:

    x(t)=x1+bxt+cxt2+dxt3y(t)=y1+byt+cyt2+dyt3

    计算机中的字母或者数字便是使用二维贝塞尔曲线画出来的,如Times Roman字体中大写T字符由下面的16条贝塞尔曲线构成,每条曲线包含x1y1x2y2x3y3x4y4来定义每一段的贝塞尔样条。

    237 620 237 620 237 120 237 120; 237 120 237 35 226 24 143 19; 143 19 143 19 143 0 143 0; 143 0 143 0 435 0 435 0; 435 0 435 0 435 19 435 19; 435 19 353 23 339 36 339 109; 339 109 339 108 339 620 339 620; 339 620 339 620 393 620 393 620; 393 620 507 620 529 602 552 492; 552 492 552 492 576 492 576 492; 576 492 576 492 570 662 570 662; 570 662 570 662 6 662 6 662; 6 662 6 662 0 492 0 492; 0 492 0 492 24 492 24 492; 24 492 48 602 71 620 183 620; 183 620 183 620 237 620 237 620;

    画出来的曲线为:

    image-20231103231902200

    python代码: BezierCurve[贝塞尔曲线].py

     

    最小二乘

    最小二乘与法线方程

    在方程组一章中,给出了当方程解存在时,如何找到Ax=b的解。我们将学习当解不存在的时候该怎么办。当方程不一致时,有可能方程的个数超过未知变量的个数,答案是找到第二可能好的解:即最小二乘近似.

    在插值一章中给出了如何找出多项式,并精确拟合数据点。但是如果有大量的数据点,或者采集的数据点具有一定误差,使用高阶多项式精确拟合一般不是个好方法。在这种情况下,使用简单模型 近似数据是一种更合理的方式。

    求解不一致的方程组以及近似拟合数据这两个问题一同驱动着最小二乘的发展.

     

    不一致的方程组

    如果一个方程组无解,它被称为不一致。方程无解意味着什么?可能系数不十分精确,在很多情况下,方程的个数比未知变量的个数要多,使得解不可能满足所有的方程。但是还可以通过找出与解最相似的向量x。

    如果我们选择“相似度”意味着欧氏距离相近,有一个直接的方法找到最接近的x,这个特殊的x称为 最小二乘解。

    对于Ax=b这个方程而言,如果b不在Ax可以确定的平面(A中的每一列看成基向量,对应的x的一行(就是一个值)便可看成该方向的值)上,那么就是无解的,不过可以在这个平面上找到与b最接近的向量来近似。

    最小二乘基于正交,从一点到一个平面的最短距离,由一个到平面的正交线段表示,法线方程可以确定该线段,这表示最小二乘的误差。

    假设这个近似解为x¯,则bAx¯代表余弦误差,让误差与平面正交,即

    AT(bAx¯)=0ATAx¯=ATb

    此时这个方程便是可解的,不过解出的结果与实际结果可能存在一定差距,有三种方式可以表示余项的大小。

    (1)欧式长度:||r||2=r12++rm2

    (2)2范数,平方误差:SE=r12++rm2

    (3)平均平方根误差:RMSE=SE/m=(r12++rm2)/m

    三种表达紧密相关。

     

    数据的拟合模型

    最小二乘数据拟合

    给定一组数据点(t1,y1),...,(tm,ym)

    step1 选择模型,确定用于拟合数据的参数模型,例如y=a+bt

    step2 将数据点代入模型,每个数据点生成一个方程,其中的未知变量是模型参数,例如线性模型中的c1c2为参数,这得到系统Ax=b,其中未知变量x表示未知参数。

    step3 参数的最小二乘解是法线方程ATAx¯=ATb的解。

    最小二乘是数据压缩的经典例子,输入包含一组数据点,输出是一个尽可能好的数据拟合模型,其中具有较少的参数。通常,使用最小二乘的原因是使用合理的底层模型替换噪声数据,该模型通常用于信号预测以及分类。

    最小二乘的条件

    ATA的条件数过大,难以在双精度算术中处理,即使原始的最小二乘问题条件还可以,但是这个法线方程是病态问题。

     

    模型概述

    周期数据

    周期数据使用周期模型。例如外层大气温度遵循大时间尺度的循环,包含由地球自转和绕太阳公转控制的每天和每年的循环。

    如拟合一天为周期的温度变化,可以使用模型y=c1+c2cos2πt+c3sin2πt+c4cos4πt

    数据线性化

    以描述人口变化的指数模型为例:

    y=c1ec2t

    不能直接进行最小二乘拟合,这是由于c2在模型方程中不是线性的。一旦将数据点代入模型,求解的困难很明显:要求解的方程系数不是线性,不能表示为线性方程组Ax=b。因而我们导出的法线方程不能使用。

    不过有两种方法可以用来解决这个问题,第一种方法是直接求解非线性最小二乘方程(在后面介绍),第二种是通过取对数来使问题线性化。

    lny=ln(c1ec2t)=lnc1+c2t

    同样对于幂法则模型y=c1ec2t,也可以采用这种方法解决。

     

    QR分解

    格拉姆-施密特正交与最小二乘

    格拉姆-施密特方法是对一组向量正交化。给定一组输入的m维向量,目的是找出正交坐标系统,获取由这些向量张成的空间。更精确地讲,给定n个线性无关的输入向量,该方法计算n个彼此垂直的单位向量(单位长度由2范数定义),张成和输入向量相同的子空间。

    格拉姆-施密特正交的思想是对向量减去其投影在其它已经正交好的向量上的分量,得到正交于这些正交好的向量,在进行单位化即可。令rij=qiTAj表示第j个向量投影在第i个已正交化好的向量上的分量长度,再乘上qi(表方向的单位向量)即表示第j个向量在第i个已正交化好的向量上的投影。定义rjj=||yj||2rij=qiTAj,则可以将格拉姆-施密特正交换一种方式表达

    A1=r11q1A2=r12q1+r22q2Aj=r1jq1++rj1,jqj1+rjjqj

    或者写成矩阵形式

    (A1An)=(q1qn)[r11r12r1nr22r2nrnn]

    或A=QR,其中A是包含列向量Aj的矩阵,我们将其称为消减QR分解。关于Aj线性无关的假设保证主对角线系数rjj非0。

    经典格拉姆-施密特正交

    Aj(j=1,...,n)为线性无关向量。

    for j = 1,2,...,n

    y=Aj

    for i = 1,2,...,j-1

    rij=qiTAj

    y=yrijqi

    end

    rjj=||y||2

    qj=y/rjj

    end

    当成功分解后,通常会填满正交单位向量的矩阵,从而得到Rm完整的基。可以通过在Aj中加入m-n个额外的向量,因而m向量可以张成Rm,并实现格拉姆-施密特方法。如下所示

    (A1An)=(q1qm)[r11r12r1nr22r2n00rnn00]

    注意A是m×n矩阵,Q是m×m方阵,上三角矩阵R是m×n矩阵,m>n,这称为A的完全QR分解,Q在数值分析中具有特殊的地位,并具有一个特殊的定义。

    定义:当QT=Q1,方阵Q正交。

    引理:如果Q是m×m正交矩阵,x是m维向量,则

    ||Qx||2=||x||2

    LU分解是对高斯消去中信息进行有效编码的方式。以相同方式,QR分解记录了矩阵正交化的信息 ,即构造一个正交集,张出由A的列向量构成的空间。使用正交矩阵计算更好的原因是

    (1)根据定义它们很容易求逆,

    (2)由引理知,它们不会放大误差。

    QR分解可用于求解n个方程n个未知量的方程组,但是计算复杂度比LU分解大。

    QR分解可以用于最小二乘,令A是m×n矩阵,其中m≥n。为了最小化||Axb||2,使用上面的引理得到:

    Qx2=x2QRxb2=Q1(QRxb)2=RxQTb2

    欧几里得范数中的向量是

    [e1enen+1em]=[r11r12r1nr22r2n00rnn00][x1xn][d1dndn+1dm]

    其中d=QTb,假设rii0,则误差向量e上面的部分(e1,,en)可以通过回代变为0,选择xi使得误差向量下面部分没有受到影响;显然(en+1,,em)=(dn+1,,dm)。因而,通过使用上面部分回代得到的x可以最小化最小二乘的解,最小二乘误差是||e||22=dn+12++dm2

    通过QR分解实现最小二乘

    给定m×n不一致系统

    Ax=b

    找出完全QR分解A=QR,令

    R^=R的上n×n子矩阵

    d^=d=QTb的上面的n个元素

    求解R^x¯=d^得到最小二乘解x¯

    相比于法线方程求解最小二乘,这个更精确。,可以更好处理病态问题。

     

    改进的格拉姆-施密特正交

    改进的格拉姆-施密特正交

    Aj(j=1,...,n)为线性无关向量。

    for j = 1,2,...,n

    y=Aj

    for i = 1,2,...,j-1

    rij=qiTy

    y=yrijqi

    end

    rjj=||y||2

    qj=y/rjj

    end

    与经典格拉姆-施密特唯一不同的是,Aj被y在最内层循环中替换。从几何上来讲,例如当减去Ajq2方向的投影时,应该减去Aj的余项y的投影,在余项中q1部分已经减掉了,而不是将Aj自己投在q2上。

     

    豪斯霍尔德反射子

    尽管改进的格拉姆-施密特正交是计算矩阵的QR分解的有效方式,但是这不是最好方式。另一种方法是使用豪斯霍尔德反射,这种方式需要更少的计算,同时在舍入误差放大的意义上来讲这种方法也更稳定。

    豪斯霍尔德反射子是正交矩阵,通过m-1维平面反射m维向量。这意味着当每个向量乘上矩阵后,长度保持不变,使得豪斯霍尔德反射被称为移动向量的完美形式。给定一个向量x,我们要重新找出一个相同长度的向量w,计算豪斯霍尔德反射得出矩阵H满足Hx=w。

    如图所示,原始方法非常清晰。画出m-1维平面二分x和w,并和连接它们的向量垂直,然后通过该平面反射所有向量。

    image-20231104104642216

    引理:假设x和w是具有相同欧几里得长度的向量,||x||2=||w||2,则w-x和w+x正交。

    定义向量v=w-x,考虑投影矩阵

    P=vvTvTv

    投影矩阵是一个满足P2=P的矩阵。并且P满足Pv=v。几何上说,对于任何向量u,Pu是u在v上的投影。从上图可以看到,如果从x中两次减去Px,就可以得到w。令H=I-2P,则

    Hx=x2Px=wv2vvTxvTv=wvvvTxvTvvvT(wv)vTv=wvvTvvTvvvTxvTvvvT(wv)vTv=wvvTwvTvvvT(wv)vTv=wvvT(2wv)vTv=wv(wx)T(w+x)vTv=w

    矩阵H被称为豪斯霍尔德反射子,注意到H是对称并正交的矩阵。

    定理(豪斯霍尔德反射子):令x和w是向量,||x||2=||w||2,并定义v=w-x,则H=I2vvT/vTv是对称正交矩阵,并且Hx=w。

    使用豪斯霍尔德反射子推导出QR分解的一个新方式,使用反射子的目的是将列向量x移动到坐标轴,并以此将0放在矩阵中。

    我们从矩阵A开始,希望写做形如A=QR的方程。令x1是A的第一列,令w=±(x12,0,...,0)为第一个坐标轴上的向量,它们的欧几里得长度相同。为了数值稳定,一般将x的第一个元素选为正号,以避免在生成v的过程中出现两个近似相等的数字相减。生成豪斯霍尔德反射子并满足H1x=w,在4×3的情况下,用A乘上H1得到

    H1A=H1[××××××××××××]=[×××0××0××0××]

    我们已经在A中引入一些0,我们希望继续这种方式直到A变为上三角;然后我们将得到QR分解中的R。对于一个3列的矩阵,还需要做两次这样的操作,即

    H3H2H1A=RA=H11H21H31R=H1H2H3R=QR

    举例:使用豪斯霍尔德反射子找到矩阵A的QR分解

    A=[142322]

    豪斯霍尔德反射子H1将第一列x=[1,2,2]移动到w=[||x||2,0,0]=[3,0,0],v=w-x=[2,-2,-2]

    H1=I2vvTvTv=[100010001]212[444444444]=[132323231323232313]
    H1A=[132323231323232313][142322]=[320304]

    下面使用H^2将向量 x^=[-3,-4]移动w^=[5,0],

    H^2=I2vvTvTv=[1001]280[64323216]=[0.60.80.80.6]H^2A^2=[0.60.80.80.6][34]=[50]

    因此,有

    H2H1A=[10000.60.800.80.6][132323231323232313][142322]=[320500]=R
    A=H1H2R=QR

    与格拉姆-施密特正交相比,QR分解的计算代价更低,可以更好地得到单位正交向量以及具有更低的内存需求。

     

    广义最小余项(GMRES)方法

     

    Krylov方法

    GMRES属于Krylov方法,这些方法依赖精确的Krylov空间的计算,该空间是向量{r,Ar,...,Akr}所张的空间,其中r=bAx0是初始估计的余项。由于向量Akr对于大的k倾向于一个共有方向,Krylov空间的基必须认真计算。找出Krylov空间精确的基需要正交化计算方法,例如格拉姆-施密特或者豪斯霍尔德反射子方法。

    GMRES的思想是在特殊矢量空间,即Krylov空间寻找初始估计的x0的改进,Krylov空间由余项r和它与非奇异矩阵A的积所张成。在该方法的第k步,加入Akr以扩大Krylov空间,重新对基进行正交化,然后通过最小二乘获取改进并加入到x0中。

    广义最小余项方法

    x0=初始估计

    r=bAx0

    q1=r/||r||2

    for k = 1,2,...,m

    y=Aqk

    for j = 1,2,...,k

    hjk=qjTy

    y=yhjkqj

    end

    hk+1,k=||y||2(如果hk+1,k=0,跳过下一行,并在底端终止)

    qk+1=y/hk+1,k

    最小化 ||Hck[||r||2,0,0,...,0]T||2得到ck

    xk=Qkck+x0

    end

    迭代的xk是系统Ax=b的近似解,如果hk+1,k=0,表示步骤k是最后一步。因为xk在循环中没有使用,最小化得到ck甚至可以放到循环外最后来做。如果条件数是一个重要问题,在内环进行的格拉姆-施密特正交步骤可以用豪斯霍尔德正交替换。GMRES的典型用途是用于大规模稀疏的n×n矩阵A。 理论上,算法经过n步终止,只要A是非奇异矩阵就可以得到解x。在大多数情况下,目标是仅仅运行k步,k比n要小得多。注意到矩阵Qk是n×k矩阵,并不保证是稀疏矩阵。因而内存也可能限制GMRES方法中的步数k。如果k步迭代后没能足够趋近解,而且如果n×k矩阵Qk变得大得难以处理,有一个简单想法:扔掉Qk重新开始GMRES方法,使用当前的最优估计xk作为新的x0,这种方法被称为重启GMRES

     

    预条件GMRES

    在预条件GMRES背后的概念和共枙梯度法非常相似,从非对称的线性方程组Ax=b开始,我们试图求解M1Ax=M1b,其中M是方程组中讨论的预条件子。

    预条件GMRES

    x0=初始估计

    r=M1(bAx0)

    q1=r/||r||2

    for k = 1,2,...,m

    w=M1Aqk

    for j = 1,2,...,k

    hjk=wTqj

    w=whjkqj

    end

    hk+1,k=||w||2

    qk+1=w/hk+1,k

    最小化 ||Hck[||r||2,0,0,...,0]T||2得到ck

    xk=Qkck+x0

    end

     

    非线性最小二乘

    高斯-牛顿方法

    首先会议多变量的牛顿方法,并将其应用在列向量函数F(x)T=(rTDr)T=(Dr)Tr

    DF(x)T=D((Dr)Tr)=(Dr)TDr+i=1mriDci

    其中ci是Dr的第i列,注意到Dci=Hriri的2阶偏导数矩阵或称为海森矩阵如下:

    Hri=[2rix1x12rix1xn2rixnx12rixnxn]

    高斯牛顿方法

    为了最小化

    r1(x)2++rm(x)2

    x0=初始向量,

    for k = 0,1,2,...

    A=Dr(xk)

    ATAvk=ATr(xk)

    xk+1=xk+vk

    end

    Dr(x,y)=[rx1ry1rxnryn]

    具有非线性参数的模型

    定义矩阵Dr是误差ri关于参数cj偏导数矩阵

    (Dr)ij=ricj=fcj(ti)

    举例:拟合模型c1ec2t

    定义余项

    r=[c1ec2t1y1c1ec2tmym]

    并相对c1c2计算偏导数得到

    r=[c1ec2t1y1c1ec2tmym]Dr=[ec2t1c1t1ec2t1ec2tmc1tmec2tm]

    最小二乘问题中的非线性带来额外的挑战。法线方程以及QR方法只要系数矩阵A满秩都可以找到唯一解。而对于非线性问题的高斯-牛顿迭代可能收敛到多个极小平方误差中的一个。尽可能使用初始向量的合理近似,有助于收敛到绝对极小。

     

    Levenberg-Marquardt方法

    非线性最小二乘法遇到病态系数矩阵时,会面临巨大的挑战,Levenberg-Marquardt方法使用“正则化项”部分修复这个问题。

    Levenberg-Marquardt方法

    最小化

    r1(x)2++rm(x)2

    x0=初始向量,λ=常数

    for k = 0,1,2,...

    A=Dr(xk)

    (ATA+λdiag(ATA))vk=ATr(xk)

    xk+1=xk+vk

    end

    提高正则化参数λ可以增强矩阵ATA对角线元素的作用,这通常可以改善条件数,允许方法从一个比高斯-牛顿更宽的初始估计x0开始,并实现收敛。

     

    数值微分和积分

    计算微积分解决的主要问题是计算函数的导数和积分。

    数值微分

    有限差分公式

    二点前向差分公式

    f(x)=f(x+h)f(x)hh2f(c)

    其中c在x和x+h之间,h2f(c)是误差项。

    定理(推广中值定理):令f是区间[a,b]上的连续函数,令x1,...,xn是区间[a,b]上的点,a1,...,an>0,则在a和b之间存在数字c满足

    (a1++an)f(c)=a1f(x1)++anf(xn)

    三点中心差分公式

    f(x)=f(x+h)f(xh)2hh26f(c)

    其中x-h<c<x+h。

    二阶导数的三点中心差分公式

    f(x)=f(xh)2f(x)+f(x+h)h2h212f(4)(c)

    其中x-h<c<x+h。

    当h->0时,误差会先减小再增大,后增大的原因是计算时有效数字的丢失。对于三点中心差分公式而言,f'(x)的及其近似误差的绝对值上界是

    E(h)h26f(c)+εmachh

    外推

    假设我们有n阶公式F(h)近似一个给定量Q,这个阶数意味着

    QF(h)+Kh

    其中K大约是我们感兴趣的h区间上的一个常数,相关的例子是

    f(x)=f(x+h)f(xh)2hf(ch)6h2

    ch的位置在x+h和x-h之间,只要f合理光滑,h不太大,那么f(ch)f(x)相差不会太大。这样我们就可以把n阶公式变换成更高阶的公式,我们期望

    QF(h/2)12n(QF(h))

    n阶公式的外推

    Q2nF(h/2)F(h)2n1

    举例:应用外推公式

    F4(x)=22F2(h/2)F2(h)221=[4f(x+h/2)f(xh/2)hf(x+h)f(xh)2h]/3=f(xh)8f(xh/2)+8f(x+h/2)f(x+h)6h

    数值积分的牛顿-科特斯公式

    梯形法则

    x0x1f(x)dx=h2(y0+y1)h312f(c)

    其中h=x1x0,c在x0x1之间。

     

    辛普森法则

    x0x2f(x)dx=h3(y0+4y1+y2)h590f(4)(c)

    其中h=x2x1=x1x0,c在x0x2之间。

    定义:数值积分方法的精度是最大的整数k,使用该积分方法可以得到所有k阶或者更低阶多项式积分的精确值。

    梯形法则的精度是1,辛普森法则的精度是3。

    复合牛顿-科特斯公式

    由于积分在区间的所有子区间上具有可加性,我们可以通过除法把整个区间变为很多小区间再计算积分,在每个小区间上使用法则,然后再求和。这种策略被称为复合数值积分

    复合梯形法则

    abf(x)dx=h2(y0+ym+2i=1m1yi)(ba)h212f(c)

    其中h=(b-a)/m,c在a和b之间。

    复合辛普森公式

    abf(x)dx=h3(y0+y2m+4i=1my2i1+2i=1m1y2i)(ba)h4180f(4)(c)

    其中c在a和b之间。

    开牛顿-科特斯方法

    中点法则

    x0x1f(x)dx=hf(w)+h324f(c)

    其中h=(x1x0),w是中点x0+h/2,c在x0x1之间。

    中点法则也可用于减少所需函数求值的次数与梯形法则相比,闭牛顿-科特斯方法具有相同的阶数,只需要一次函数求值,而不是两次。而且误差项是梯形法则误差项大小的一半。

    复合中点法则

    x0x1f(x)dx=hi=1mf(wi)+(ba)h224f(c)

    其中h=(b-a)/m,c在a和b之间,wi是[a,b]中m个相等子区间的中点。

     

    龙贝格积分

    龙贝格积分是对复合梯形法则应用外推的结果,给定近似M的法则N(h),该法则依赖步长h,如果知道法则的阶,则可以对法则进行外推。

    龙贝格积分

    R11=(ba)f(a)+f(b)2

    for j = 2,3,...n

    hj=ba2j1

    Rj1=12Rj1,1+hji=12j2f(a+(2i1)hj)

    for k = 2,...,j

    Rjk=4k1Rj,k1Rj1,k14k11

    end

    end

    最后的结果便是Rjj,龙贝格积分的常用停止条件是计算新的一行直到相邻的对角线元素Rjj差异小于当前的容差。

     

    自适应积分

    到目前我们学到的积分近似方法都使用相等的步长。一般来说,小步长可以提高精度变化剧烈的函数将需要更多步,因而带来更多的计算时间,这是由于需要更小的步长跟踪函数所有的变化。函数可能在定义域某些部分变化剧烈,而在另外一些部分变化缓慢,满足前一部分的容差的步长在后一部分就显得过于精细了。

    通过使用积分误差公式,可以在运算中推出一个标准,其步长对于特定的子空间适合。这种方法称为自适应积分

    自适应积分

    对给定容差TOL近似abf(x)dx

    c = (a+b)/2

    S[a,b]=(b-a)(f(a)+f(b))/2

    if S[a,b]S[a,c]S[c,b]<3TOL(baborigaorig)

    接受S[a,c]+S[c,b]作为区间[a,b]上的近似

    else

    对于区间[a,c]和[c,b]递归重复上面的步骤

    end

    其中tol0是自己设置的参数,如设置为0.005。

     

    高斯积分

    定义:一组在区间[a,b]上的非0函数{p0,...,pn}在区间[a,b]正交,当

    abpj(x)pk(x)dx={0jk0j=k

    定理:如果{p0,...,pn}是区间[a,b]上多项式的正交集,其中degpi=i,则{p0,...,pn}是在区间[a,b]由区间[a,b]由至多n个多项式所张的向量空间的基。

    定理:如果{p0,...,pn}是区间[a,b]上多项式的正交集,且degpi=i,则pi在区间(a,b)上有i个不同的根。

    高斯积分

    11f(x)dxi=1ncif(xi)

    其中

    ci=11Li(x)dx,i=1,,n
    Li=(xx1)(xxi)(xxn)(xix1)(xixi)(xixn)

    上面带横线的项表示不进行计算,ci可以查表获得(包括n=1,2,3,4)

    定理:高斯积分方法,在区间[—1,1]上使用n阶勒让德多项式,具有2n-1阶精度。

    一般区间[a,b]的积分可以通过替代t=(2x-a-b)/(b-a)得到

    abf(x)dx=11f((ba)t+b+a2)ba2dt

    常微分方程

    初值问题

    欧拉方法

    通过跟随箭头计算“求解”微分方程。从初始条件(t0y0)开始,然后沿着在那里指定的方向。在移动很小的距离后,在新点(t1y1)重新计算斜率,根据新的斜率继续移动,重复这个过程。

    欧拉方法

    w0=y0

    wi+1=wi+hf(ti,wi)

    h是步长,如取0.2,y'=f(t,w)。

    举例:用欧拉方法求解如下问题

    {y=ty+t3y(0)=y0t[0,1]

    欧拉方法则对应如下迭代,h设置为0.1

    w0=1wi+1=wi+h(tiwi+ti3)

    解的存在性、唯一性和连续性

    定义:当存在常数L(称为利普希茨常数)对矩形 S=[a,b]×[α,β]中的每对(t,y1),(t,y2)满足

    |f(t,y1)f(t,y2)|L|y1y2|

    函数f(t,y)相对于变量y在S上利普西茨连续

    定理:假设f(t,y)定义在集合[a,b]×[α,β]并且α<ya<β,函数对于变量y是利普希茨连续,则在a与 b之间存在c,使得初值问题

    {y=f(t,y)y(a)=yat[a,c]

    有唯一解y(t)。而且,如果f在[a,b]×(-∞,∞)是利普西茨连续,则它在[a,b]上存在唯一解。

    定理:假设f(t,y)在集合S=[a,b]×[α,β]上关于y是连续的,如果Y(t)和Z(t)是微分方程

    y=f(t,y)

    在S上的解,分别具有初值条件Y(z)和Z(a),则

    |Y(t)Z(t)|eL(ta)|Y(a)Z(a)|

    一阶线性方程

    一组容易求解的特定的常微分方程提供了ODE求解问题中具有指导意义的例子,这组方程是一阶方程,其右侧是关于y的线性函数。考虑初值问题

    {y=g(t)y+h(t)y(a)=yat[a,b]

    首先注意到如果g(t)在区间[a,b]上连续,唯一解存在,使用L=max[a,b]g(t)作为利普希茨常数可以使用技巧找到解,即对方程乘上“积分因子”。

    积分因子是eg(t)dt,在方程两侧同时乘这个因子得到

    (yg(t)y)eg(t)dt=eg(t)dth(t)(yeg(t)dt)=eg(t)dth(t)yeg(t)dt=eg(t)dth(t)dty(t)=eg(t)dteg(t)dth(t)dt

    IVP求解器的分析

    显式梯形方法

    显式梯形方法

    w0=y0

    wi+1=wi+h2(f(ti,wi)+f(ti+h,wi+hf(ti,wi)))

    泰勒方法

    泰勒方法基本思想是直接利用泰勒展开,假设解y(t)是(k+1)阶连续可微的函数。给定在解曲线上的当前点(t,y(t)),目标是对于某个步长h,使用微分方程的信息,用y(t)来表达y(t+h)。

    K阶泰勒方法

    w0=y0

    wi+1=wi+hf(ti,wi)+h22f(ti,wi)++hkk!f(k1)(ti,wi)

    举例:对于一阶线性方程确定二阶泰勒方法

    {y=ty+t3y(0)=y0

    由于f(t,y)=ty+t3

    f(t,y)=ft+fyf=y+3t2+t(ty+t3)

    所以二阶泰勒方法的迭代公式为

    wi+1=wi+h(tiwi+ti3)+h22(wi+3ti2+ti(tiwi+ti3))

    常微分方程组

    高阶方程

    单个的高阶微分方程可以转化为一个方程组,令

    y(n)=f(t,y,y,y,,y(n1))

    是n阶常微分方程,定义新的变量

    y1=yy2=yy3=yyn=y(n1)

    注意到原始的常微分方程写成

    yn=f(t,y1,y2,,yn)

    两者放在一起

    y1=y2y2=y3y3=y4yn1=ynyn=f(t,y1,y2,,yn)

    这些方程将n阶的微分方程转化为一阶微分方程组,该方程组可以使用欧拉或者梯形方法求解。

     

    龙格-库塔方法和应用

    龙格-库塔家族

    我们已经知道一阶的欧拉方法以及二阶的梯形方法。除了梯形方法外,还有其他龙格-库塔类型的二阶方法。

    中点方法

    w0=y0

    wi+1=wi+hf(ti+h2,wi+h2f(ti,wi))

    4阶龙格-库塔方法(RK4)

    wi+1=wi+h6(s1+2s2+2s3+s4)

    其中

    s1=f(ti,wi)s2=f(ti+h2,wi+h2s1)s3=f(ti+h2,wi+h2s2)s4=f(ti+h,wi+hs3)

     

    可变步长方法

    直到现在,步长h在ODE求解器实现中一直是一个常数。但是,并没有原因说h在求解过程中不能改变。我们希望改变步长的一个原因是,问题的解会在缓慢变化的周期和迅速变化的周期之间移动。将固定步长变得足够小使其精确跟踪快速变化,意味着解的其他部分求解都会是令人难以忍受的缓慢。

    龙格-库塔嵌入对

    可变步长方法的关键思想是检测当前步生成的误差。用户设置的容差在当前步必须能够满足.然后设计方法(1)如果超出容差,必须拒绝误差并将步长减小,或者(2)如果满足容差,接受步长并选择对于下一步适合的步长h。关键是近似在每步中产生的误差。首先假设我们已经找到这样的方式并解释如何改变步长。

    有RK2/3(龙格-库塔2阶/3阶嵌入对)、RK4/5,不过多赘述。

     

    隐式方法和刚性方程

    后向欧拉方法

    w0=y0

    wi+1=wi+hf(ti+1,wi+1)

    举例:对初值问题使用后向欧拉方法

    {y=y+8y29y3y(0)=1/2t[0,3]

    使用后向欧拉方法

    wi+1=wi+hf(ti+1,wi+1)=wi+h(wi+1+8wi+129wi+13)

    z=wi+1,必须求解方程z=wi+h(z+8z29z3),或者对于未知的z求解

    9hz38hz2+(1h)zwi=0

    可以使用牛顿方法进行估计

    znew=z9hz38hz2+(1h)zwi27hz216hz+1h

    znew替换z并重复进行,对于每个后向欧拉步,使用牛顿方法直到znewz小于容差(比近似微分方程解的误差还要小)。

     

    多步方法

    我们已经研究过的龙格-库塔方法家族包含单步方法,意味着新一步的wi+1可基于微分方程和上一步的wi得到。 多步方法则给出不同的方式:使用除了wi之外的知识帮助下一步的求解。这将带来OD E求解器的阶数和单步方法的阶数相同,但是大量必要的计算会被替换为求解过程中已经计算的值 的插值。

    构造多步方法

    Adams-Bashforth两步方法

    wi+1=wi+h[32f(ti,wi)12f(ti1,wi1)]

    显式多步方法

    Adams-Bashforth三步方法(三阶)

    wi+1=wi+h12[23fi16fi1+5fi2]

    Adams-Bashforth四步方法(四阶)

    wi+1=wi+h24[55fi59fi1+37fi29fi3]

    隐式多步方法

    隐式梯形方法(二阶)

    wi+1=wi+h2[fi+1+fi]

    Adams-Moulton两步方法(三阶)

    wi+1=wi+h12[5fi+1+8fifi1]

    Milne-Simpson方法

    wi+1=wi+h3[fi+1+4fi+fi1]

    Adams-Moulton三步方法(四阶)

    wi+1=wi+h24[9fi+1+19fi5fi1+fi2]

    Adams-Moulton四步方法(五阶)

    wi+1=wi+h720[251fi+1+646fi246fi1+106fi219fi3]

     

    边值问题

    有限差分方法

    线性边值问题

    令y(t)是一个至少四阶连续的函数,可以得到一阶导数的离散近似为

    y(t)=y(t+h)y(th)2hh26y(c)

    二阶导数的近似为

    y(t)=y(t+h)y(th)2hh26y(c)y(t)=y(t+h)2y(t)+y(th)h2+h212y(c)

    有限差分方法包含使用离散版本替换微分方程中的导数,并求解得到的更简单的代数方程获取正确值yi的近似wi,边界条件则在需要的时候代入方程组中。

    举例:使用有限差分求解

    {y=4yy(0)=1y(1)=3

    使用二阶导数的中心差分形式,在ti处的有限差分形式如下

    wi+12wi+wi1h24wi=0

    或者等价地

    wi1+(4h22)wi+wi+1=0

    当n=3时,区间大小h=1/(n+1)=1/4,再插入边界条件w0=1w4=3,给出方程

    {1+(4h22)w1+w2=0w1+(4h22)w2+w3=0w2+(4h22)w3+3=0

    代入h得到三对角线矩阵方程

    [941019410194][w1w2w3]=[103]

    非线性边值问题

    对于非线性边值问题,使用多变量牛顿方法进行迭代,迭代公式为

    wk+1=wkDF(wk)1F(wk)

    一般地,完成迭代的最好方式是求解方程DF(wk)Δw=F(wk)中的Δw=wk+1wk

    举例:求解非线性边值问题

    {y=yy2y(0)=1y(1)=4

    离散化后方程为

    wi1(2+h2)wi+h2wi2+wi+1=0

    其中w0=ya=1wn+1=yb=4

    函数F(w)如下

    F[w1w2wn1wn]=[ya(2+h2)w1+h2w12+w2w1(2+h2)w2+h2w22+w3wn2(2+h2)wn1+h2wn12+wnwn1(2+h2)wn+h2wn2+yb]

    排列与有限元方法

    和有限差分方法相似,排列和有限元方法背后的思想是将边值问题消减为一组可求解的代数方程。但是,并不是通过使用有限差分替换微分方程中的导数进行离散化,而是对于解给出函数形式,其对应的参数使用该方法拟合。

    选择一组基函数ϕ1(t),...,ϕn(t),它们可能是多项式、三角函数、样条,或者其他简单函数,然后考虑可能的解

    y(t)=c1ϕ1(t)++cnϕn(t)

    找出近似解的问题被简化为确定ci的值。

     

    排列

    考虑BVP问题

    {y=f(t,y,y)y(a)=yay(b)=yb

    选择n个点,从边值点a开始,在b点结束,即

    a=t1<t2<<tn=b

    排列方法将的候选解代入微分方程中,并估计微分方程在n个点上的值,得到关于n个未知变量c1,...,cn的方程,ti的值为ti=a+i1n1(ba)

    选择基函数ϕj(t)=tj1,其中1≤j≤n。解的形式如下

    y(t)=j=1ncjϕj(t)=j=1ncjtj1

    所以可以写出n个未知变量c1,...,cn的n个方程,其中第一个和最后一个方程是边值条件

    j=1ncjaj1=y(a)j=1ncjbj1=y(b)

    余下的n-2个方程来自微分方程在ti上的求值,其中2≤i≤n-1。将微分方程y"=f(t,y,y')应用到y(t)=j=1ncjtj1得到

    j=1n(j1)(j2)cjtj3=f(t,j=1ncjtj1,j=1n(j1)cjtj2)

    ti点求值得到n个方程,然后求解ci。如果微分方程是线性的,则关于ci的方程组也是线性的,并可以求解。

    在排列方法中选择三角函数作为基带来了傅里叶分析以及谱方法,这在边值方程和偏微分方程中都大量应用。这是一个“全局的方法,其中对于t的一个非常大的范围中,基函数非0,但是具有很好的正交性质。

     

    有限元以及Galerkin方法

    选择样条函数作为基函数带来有限元方法.在这种方法中,每个基函数仅在t的一个小区间上非0。有限元方法大量用于高维情况下的BVP和PDE,特别是当存在不规则的边界条件的情况下,其中使用标准的基函数做参数化变得困难。Galerkin方法最小化微分方程在解上的平方误差。这带来关于ci的一个不同的方程组。

    考虑一个BVP问题

    {y=f(t,y,y)y(a)=yay(b)=yb

    解为y,余项r=y"-f,需要使得微分方程两侧的差异尽可能小。和最小二乘方法类似,这可以通过选择y使得余项与潜在解的向量空间正交得到。

    对于区间[a,b],定义平方可积分的函数对应的向量空间

    L2[a,b]={[a,b]y(t)|aby(t)2dt}

    L2函数空间具有内积

    y1,y2=aby1(t)y2(t)dt

    <y1,y2>=0时,两个函数y1y2L2[a,b]正交。由于L2[a,b]是无限维的向量空间,通过有限计算,不能得到余项r与所有L2[a,b]正交,但是可以使用已有的计算资源,选择基和L2张得尽可能一致,令n+2个基函数的集合表示为ϕ0(t),...,ϕn+1(t)

    Galerkin方法包含两个主要思想,第一个是最小化r,其中强制其在L2内积的意义上与基函数正交,这意味着强制ab(yf)ϕi(t)dt=0,或

    abyϕi(t)dt=abf(t,y,y)ϕi(t)dt

    其中0≤i≤n+1,该形式被称为边值问题的弱形式。

    Galerkin方法的第二个思想是部分使用积分消去二阶导数,注意到

    abyϕi(t)dt=ϕi(t)y(t)|ababy(t)ϕi(t)dt=ϕi(b)y(b)ϕi(a)y(a)aby(t)ϕi(t)dt

    可以得到

    abf(t,y,y)ϕi(t)dt=ϕi(b)y(b)ϕi(a)y(a)aby(t)ϕi(t)dt

    以如下函数形式求解ci

    y(t)=i=0n+1ciϕi(t)

    Galerkin方法的两个思想使其可以方便地使用极简单的函数作为有限元ϕi(t),我们将仅仅介绍分段线性B样条函数,从t轴上的数据格点t0<t1<<tn1<tn开始,对于i=1,...,n定义

    ϕi(t)={tti1titi1ti1<ttiti+1tti+1titi<tti+10otherwise

    同时定义

    ϕ0(t)={t1tt1t0t0t<t10otherwiseϕn+1(t)={ttntn+1tntn<ttn+10otherwise

    对于ϕi(tj)而言,只有当i=j时,才为1,否则为0。对于一组数据点(ti,ci),定义分段线性B样条

    S(t)=i=0n+1ciϕi(t)

    同时有S(t)=i=0n+1ciϕi(tj)=cj

    第一个和最后一个ci通过排列方法得到

    y(a)=i=0n+1ciϕi(a)=c0ϕ0(a)=c0y(b)=i=0n+1ciϕi(b)=cn+1ϕn+1(b)=cn+1

    对于i=1,...,n,使用有限元方程

    abf(t,y,y)ϕi(t)dt+aby(t)ϕi(t)dt=0

    代入y(t)=ciϕi(t),利用ϕ函数的性质可以求解。

     

    偏微分方程

    前向差分方法

    定义偏微分方程(PDE)如下

    {ut=Duxxaxb,t0u(x,0)=f(x)axbu(a,t)=f(t)t0u(b,t)=r(t)t0

    使用数值微分和积分的离散公式近似在x和t方向的导数。例如,使用中心差分公式,关于x的二阶导数如下

    uxx(x,t)1h2(u(x+h,t)2u(x,t)+u(xh,t))

    误差为h2uxxxx(c1,t)/12,对于时间变量的一阶导数的前向差分公式得到

    ut(x,t)=1k(u(x,t+k)u(x,t))

    其中误差kutt(x,c2)/2xh<c1<x+ht<c2<t+h。带入点(xi,tj)处的偏微分方程中得到

    Dh2(wi+1,j2wij+wi1,j)1k(wi,j+1wij)

    通过时间步进可以得到(定义σ=DKh2

    wi,j+1=wij+DKh2(wi+1,j2wij+wi1,j)=σwi+1,j+(12σ)wij+σwi1,j

    写成矩阵形式

    [w1.j+1wm,j+1]=[12σσ00σ12σσ00σ12σ0σ00σ12σ][w1.jwm,j]+σ[w0,j00wm+1,j]

    后向差分方法

    此时离散的差分公式为(定义σ=DKh2

    1k(wijwi,j1)=Dh2(wi+1,j2wij+wi1,j)σwi+1,j+(1+2σ)wijσwi1,j=wi,j1

    矩阵形式为

    [1+2σσ00σ1+2σσ00σ1+2σ0σ00σ1+2σ][w1.jwm,j]=[w1.j1wm,j1]+σ[w0,j00wm+1,j]

    Crank-Nicolson方法

    Crank-Nicolson对于时间导数使用后向差分公式,并均匀组合前向差分近似和后向差分近似。

    使用如下的后向差分公式替换ut

    1k(wijwi,j1)

    使用混合差分替换uxx

    12(wi+1,j2wij+wi1,jh2)+12(wi+1,j12wi,j1+wi1,j1h2)

    可以重新组织热方程近似

    2wij2wi,j1=σ[wi+1,j2wij+wi1,j+wi+1,j12wi,j1+wi1,j1]

    或者

    σwi1,j+(2+2σ)wijσwi+1,j=σwi1,j1+(22σ)wi,j1+σwi+1,j1

     

    本章中有各种偏微分方程(抛物线方程,双曲线方程,椭圆方程,非线性偏微分方程)的解法。

     

    随机数和应用

    随机数

    伪随机数

    定义:线性同余生成器(LCG)表示为以下形式(ui为生成的[0,1]随机数)

    xi=axi1+b(modm)ui=xim

    其中a为乘子,b为偏移,m为模数。

    最小标准随机数生成器

    xi=axi1(modm)

    ui=xim

    其中m=2311a=75=16807b=0

    按照现在计算机每秒指令数的增加,这个数目显得太少了。

    randu生成器

    xi=axi1(modm)

    ui=xim

    其中a=65539=216+3m=231

    不满足随机数的独立性要求。

    最新版的matlab已经不再使用LCG作为随机数生成器,从matlab5开始,使用的是延迟斐波那契生成器。

     

    指数和正态随机数

    指数随机变量的取值满足概率分布函数p(x)=aeax,a>0,使用均匀分布随机数,可以很容易地生成指数随机数,累计分布函数为

    P(x)=Prob(Vx)=0xp(x)dx=1eax

    算法的核心思想是要找到指数随机变量的值,使得概率Prob(V≤x)是[0,1]区间上的均匀分布,即给定一个满足均匀分布的变量u之后,要使得

    u=Prob(Vx)=1eax

    求解x得到

    x=ln(1u)a

    把输入的均匀随机数u转换为指数随机数。

    对于一般的概率分布,也可以用这种方法处理。假设P(x)是需要生成的随机变量对应的概率分布函数。Q(x)=P1(x)是对应的反函数。若U[0,1]是[0,1]区间上的均匀分布的随机数,则Q(U[0,1])就是满足分布P的随机数,剩下的问题就是如何高效地计算函数Q。

    标准正态分布也可以采用反函数的方式计算,不过还有一种更有效的方法,可以同时生成两个正态分布随机数。二维标准正态分布的概率密度函数为p(x,y)=(1/2π)e(x2+y2)/2,在极坐标系下又可写成p(r)=(1/2π)er2/2的形式。由于p(r)是对极坐标轴对称的,我们只需要考虑半径r,而角度θ可以用[0,2π]区间上的均匀分布。由于p(r)是对于变量r2的指数分布,其中参数a=1/2,可以得到r2

    r2=ln(1u1)1/2

    其中u1是均匀分布随机数,进而生成两个满足标准正态分布且完全独立的随机数

    n1=rcos2πu2=2ln(1u1)cos2πu2n2=rsin2πu2=2ln(1u1)sin2πu2

    这种算法称为Box-Muller方法,此外还有一种更加高效的方法

    n1=x12ln(u1)u1n2=x22ln(u1)u1

    该方法可以避免使用三角函数,但存在一定的拒绝率(21%左右),如果x12+x22<1就接受这个值,反之重新选择。

     

    蒙特卡洛模拟

    幂律和蒙特卡洛模拟

    Ⅰ型蒙特卡罗问题是用随机采样点计算方程的均值,再乘以积分区间的体积。计算函数均值可以看成是计算一个符合同样分布函数的概率的均值。我们用E(X)表示随机变量X的期望随机变量X的方差为E[(EE(X))2],而X标准差为方差的平方根。在测量过程中,误差的期望会随着n的增加而下降,并且符合如下的规律:

    Ⅰ型或Ⅱ型蒙特卡罗模拟具有伪随机数字。

    Errorn12

    拟随机数

    拟随机数的理念是在条件允许的情况下,放弃对随机数独立性的要求。放弃独立性意味着拟随机数不再是随机的,而且也不是像伪随机数那样看上去是随机的。通过这种方式,可以使得蒙特卡罗方法能很快地收敛到正确解。拟随机数的设计是要使得生成的序列是自避(Self-avoiding)的,而不一定要保证独立性。自避定义为,在生成随机数序列的过程中,新的数会填充到比较稀疏的区域,而不会聚集到一起。

    Halton序列:设p为一个素数,例如p=2。按p进制表示写出1~n的自然数,假设第i个数字的p进制表示为bkbk1b2b1,则我们要生成的随机数就是0.b1b2bk1bk。换句话说,整个生成过程就是先写出第i个整数,然后把位数反过来写,再放到小数点的另一侧,最后输出[0,1]区间上均匀分布的第i个随机数。

    如p=3生成的前8个随机数

    i(i)3(ui)3uii(i)3(ui)3ui
    110.10.3¯5120.210.7¯
    220.20.6¯6200.020.2¯
    3100.010.1¯7210.120.5¯
    4110.110.4¯8220.220.8¯

    拟随机数的优势在于它的收敛速度更快。如果把估计误差写成计算次数n的函数的话,使用拟随机数会以更高阶的速度收敛。以下是采用拟随机数的误差收敛函数,d表示将要生成的随机数的维度:

    1型蒙特卡洛问题(拟随机数)

    Error(lnn)dn1

    2型蒙特卡洛问题(拟随机数)

    Errorn1212d

    离散和连续布朗运动

    随机游走

    随机游走模型,又称离散布朗运动。随机游走Wt定义在实数轴上,起始位置W0=0,在每个整数时刻i,都会移动一个距离si,其中si是独立同分布的随机变量。这里,我们假设每个si只能是+1或-1,且概率均为1/2。离散布朗运动定位为按如下公式生成的随机游走序列

    Wt=W0+s1+s2+cdots+st

    其中t=0,1,2,...。

    每个时刻t,Wt是一个随机变量,把这些随机变量连起来得到的{W0W1W2,...}就称为随机过程,Wt=s0+s1++st的期望为0,方差为t。

    很多随机游走模型的应用都是关于逃逸时间,又称为首次到达时间。设a,b为正整数,一个初始位置为0的随机游走序列,首次达到[-b,a]区间边缘的时刻,就称为逃逸时间。理论证明在a处(而不是-b处)逃逸的概率为b/(a+b)。对于从区间[-b,a]逃逸所需的时间,其期望值为ab。

    连续布朗运动

    标准随机游走过程在时刻t的期望为0,方差为t。假设在每个单位时间内都可以进行两次位移,且每次位移的时间为1/2个单位时间,则在时刻t的期望仍然会是0,但是方差会变为

    V(Wt)=V(s1++s2t)=2t

    这是因为一共进行了2t次位移操作。如果我们把移动次数增加到原来的k倍,我们就必须把每次移动的距离缩小到原来的1/k,这样才能保证方差不变。

    因此,Wtk表示的随机游走过程,其在水平方向上的步长为原来的1/k,在垂直方向上的步长为±1/k的等概率分布,则期望仍然为0,方差为

    V(Wtk)=i=1ktV(sik)=i=1kt[12(1k)2+12(1k)2]=kt1k=t

    当k趋向∞,极限Wt就成为了连续布朗运动,Bt=Wt是在t≥0上的随机变量,有三个重要特性

    (1)对任意t,随机变量Bt的均值为0,方差为t;

    (2)对任意t1<t2,随机变量Bt2Bt1是正态分布,且独立于Bt1的取值,实际上独立于所有Bs,0≤s≤t1

    (3)布朗运动Bt是一条连续路径。

     

    随机微分方程

    有噪声的微分方程

    定义:以实数t≥0为索引的随机变量集合xt称为时域连续随机过程

    考虑如下的有噪声的微分方程(SDE)问题:

    {dy=rdt+σdBty(0)=0

    其中r和σ为常数,这个概率过程为y(t)=rt+σBt

    SDE是以微分的形式给出的,而在ODE中使用的是导数形式。这是因为像包括布朗运动在内的很多随机过程,它们是连续但不可微的。因此,SDE

    dy=f(t,y)dt+g(t,y)dBt

    等价于积分形式的表示

    y(t)=y(0)+0tf(s,y)ds+0tg(s,y)dBs

    公式中第二个积分项称为Ito积分。

    a=t0<t1<<tn1<tn=b是在区间[a,b]上的网格点,黎曼积分定义为如下的极限

    abf(t)dt=limΔt0i=1nf(ti)Δti

    其中Δti=titi1,ti1titi。类似地,Ito积分定义为如下的极限:

    abf(t)dBt=limΔt0i=1nf(ti1)ΔBi

    其中ΔBi=BtiBti1,是布朗运动中单次位移的距离。

    布朗运动Bt的微分dBt称为白噪声。

    Ito公式

    若y=f(t,x),则

    dy=ft(t,x)dt+fx(t,x)dx+122fx2(t,x)dxdx

    有许多数值方法用于求解SDE,有Euier-Maruyania方法,Milstein方法,一阶随机龙格-库塔方法,不过多赘述。

     

     

    三角插值和FFT

    DFT和快速傅里叶变换的内容不过多赘述。

    三角插值

    DFT插值定理

    设[c,d]为参数区间,n为正整数。定义Δt=(dc)/ntj=c+jΔt,其中j=0,...,n-1。这是在区间[c,d]上均匀分布的采样点,向量x为傅里叶变换的输入,xj可以理解为一个信号的第j 个分量。例如,我们可以把x想象成为一系列的测量,每个xj都测量一个离散的采样点tj

    设y=Fnx是x的傅里叶变换,由于x是y的逆傅里叶变换,x可以表示为

    xj=1nk=0n1yk(wk)j=1nk=0n1ykei2πkj/n=k=0n1ykei2πkj/nn=k=0n1ykei2πk(tjc)dcn

    以上公式可以看成是对采样点(tj,xj)的插值,其中使用的是三角函数作为基函数,而插值系数为yk。而yk=Fnx,所以可以说傅里叶变换Fn可以将数据{xj}变换为插值系数。

    假设yk=ak+ibk,可以将插值函数写成

    Q(t)=1nk=0n1(ak+ibk)(cos2πk(tc)dc+isin2πk(tc)dc)

    把插值函数Q(t)=P(t)+I(t)分成实部和虚部两个部分,由于xj是实数,只需用到实数部分就可以插值得到xj,实部为

    P(t)=Pn(t)=1nk=0n1(akcos2πk(tc)dcbksin2πk(tc)dc)

    下标n表示三角插值公式中的项数,有时可以把Pn称为n阶三角函数。

    当n为偶数时,有

    Pn(t)=a0n+2nk=1n/21(akcos2πk(tc)dcbksin2πk(tc)dc)+an/2ncosnπ(tc)dc

    满足Pn(tj)=xj,j=0,...,n-1。

    当n为奇数时,有

    Pn(t)=a0n+2nk=1(n1)/2(akcos2πk(tc)dcbksin2πk(tc)dc)

    FFT和信号处理

    使用三角函数计算最小二乘近似问题

    正交性和插值

    因为Fn1=Fn¯T=Fn¯F¯表示共轭),即Fn是酉阵,可以给出一种特殊形式的正交阵,可以用来进行插值计算。

    定理:令f0(t),...,fn1(t)为关于时间t的函数,t0,...,tn1为实数。设n×n矩阵

    A=[f0(t0)f0(t1)f0(tn1)f1(t0)f1(t1)f1(tn1)fn1(t0)fn1(t1)fn1(tn1)]

    为实数正交阵,对y=Ax,以下函数

    F(t)=k=0n1ykfk(t)

    是经过采样点(t0,x0),...,(tn1,xn1)的插值函数,即F(tj)=xj,j=0,...,n-1。

    令y=Ax,可以直接给出插值函数

    F(t)=1ny0+2ny1cos2π(tc)dc+2ny2sin2π(tc)dc+2ny3cos4π(tc)dc+2ny4sin4π(tc)dc++1nyn1cosnπ(tc)dc

    用三角函数进行最小二乘拟合

    DFT方法可以对区间[0,1]上的均匀采样点进行三角插值,假设n为偶数

    Pn(t)=a0n+2nk=1n/21(akcos2πktbksin2πkt)+an/2ncosnπt

    公式中的项数为n,等于采样点的数目。采样点越多,需要用到的sin和cos函数项就越多。

    当样本点数目 n 很大时,一般不会去精确地拟合模型函数. 实际上,在一般应用中,需要丢失一些信息(有损压缩)而使得模型比较简单;此外,由于采样点本身就是不准确的,因此强制要求插值函数经过所有点是不恰当的。对于这两种情况,可以使用上式进行最小二乘拟合,由于在这个模型中akbk都是线性的,所以可以用法线方程求解。

    令n表示样本点xj的数目,xj是[0,1]区间上的均匀时刻tj=j/n的采样点。我们要引入正偶数m, 表示在最小二乘拟合中要用到的基函数个数,即要拟合到前m个基函数,f0(t),...,fm1(t),拟合的结果如下式

    Pm(t)=k=0m1ckfk(t)

    当m<n时,就变成了一个压缩问题,我们希望用Pm在最小平方误差的标准下拟合采样点。

    最小二乘问题需要找到系数c0,...,cm1,使得方程

    k=0m1ckfk(tj)=xj

    的误差尽可能小,用矩阵形式表示,为

    AmTc=x

    其中Am是矩阵A的前m行,假设AmT的各列是相互正交的,对于变量c的法线方程

    AmAmTc=Amx

    AmAmT是单位阵。因此,最小二乘的解为

    c=Amx

    n个点的插值多项式可以表示为

    Fn(t)=k=0n1ykfk(t)

    则解可以表示为

    Fm(t)=k=0m1ykfk(t)

    这个结论说明,对给定的n个数据点,用m<n个三角函数进行最小二乘拟合,只需计算n阶插值,再取前m个项即可。即x的插值系数Ax,在丢弃高频分量后,仍然是x的最佳拟合函数。n阶展开里的前m项,就是使用m个低阶频率可以达到的最佳拟合。这个性质反映了基函数的“正交性”。

    推论:令[c,d]为一区间,m<n为正偶数,x=(x0,...,xn1)为n维实数向量,tj=c+j(dc)/n,j=0,...,n-1。令{a0a1b1a2b2,...,an/21bn/21an/2}=Fnx为x的插值系数,且

    Pn(t)=a0n+2nk=1n/21(akcos2πk(tc)dcbksin2πk(tc)dc)+an/2ncosnπ(tc)dc

    其中j=0,...,n-1,则有(注意Pmam/2的那一项中有2,而Pn中这一项没有2,这与原书423页的引理10.10有关)

    Pm(t)=a0n+2nk=1m/21(akcos2πk(tc)dcbksin2πk(tc)dc)+2am/2ncosnπ(tc)dc

    是数据(tj,xj),j=0,...,n-1的m阶最优最小二乘拟合。

     

    声音、噪声和滤波

    滤波的用法有两种。第一种是用另一个简单的函数尽可能地拟合原始音频.这是信号的压缩,我们可以只储存m个低频分量。滤波的另一个重要应用是去噪,在一个音乐文件里,音乐和语音中可能会混有高频噪声,需要消除这些高频分量来提高音频的质量。

     

    维纳滤波

    设c为一个没有噪声的音频信号,与一个同样长度的向量r相加。得到的结果x=c+r是否是有噪声的 ?若r=c,我们认为r不是噪声,因为相加的结果只是增大c的音量,但信号是同样干净的。根据定 义,噪声是与信号无关的,即若r为噪声,则内积cTr的期望应为0。

    在一个典型应用中,我们需要从含噪声的信号x中恢复出c。信号c可能是一个重要的系统变量,但是在噪声环境下进行测量。或者如下面例子,c可能是音频采样,我们想从中去除噪声。在20世纪中期Norbert Wiener建议,以最小平方误差为目标,寻找一种最优滤波器来去掉x中的噪声。他建议寻找一个对角阵Φ,使得公式

    F1ΦFxc

    的模尽可能小,其中F表示DFT变换。它的核心思想是通过傅里叶变换后,在频域进行乘以Φ,再进行逆傅里叶变换。这又称为频域滤波,因为我们处理的是经过傅里叶变换的信号x而不是x本身。为了找到最佳的对角阵Φ,注意到

    F1ΦFxc2=ΦFxFc2=ΦF(c+r)Fc2=(ΦI)C+ΦR2

    令C=Fc和R=Fr为傅里叶变换,同时注意到根据噪声的定义

    C¯TR=FcTFr=cTF¯TFr=cTr=0

    因此模可以写成

    (ΦI)C+ΦR2=((ΦI)C+ΦR)T((ΦI)C+ΦR)=(C¯T(ΦI)+R¯TΦ)((ΦI)C+ΦR)C¯T(ΦI)2C+R¯TΦ2R=i=1n(ϕi1)2|Ci|2+ϕi2|Ri|2

    为了确定对角元素ϕi能使上式最小,对每个ϕi求导,得到

    2(ϕi1)|Ci|2+2ϕi|Ri|2=0

    对每个i,解得ϕi

    ϕi=|Ci|2|Ci|2+|Ri|2

    但一般情况下我们不知道C或R,所以需要某种近似方法,令X=Fx为傅里叶变换,沿用关于信号与噪声的独立性假设,近似地有

    |Xi|2|Ci|2+|Ri|2

    可以把最优参数选为

    ϕi|Xi|2|Ri|2|Xi|2

    这里需要用到对噪声等级的先验知识,例如,如果噪声是无关的高斯噪声(独立于信号的正态高斯分布随机数,与信号值相加),可以把|Ri|2项替换为常数(pσ)2,其中σ为噪声的标准差,p为可调节的接近1的参数。

     

    压缩

    离散余弦变换

    一维DCT

    设n为正整数,一维的n阶DCT变换定义为n×n矩阵C,其元素为

    Cij=2naicosi(2j+1)π2n

    其中i,j=0,...,n-1,

    ai={1/2ifi=01ifi=1,,n1

    定义:设C为如下矩阵,向量x=[x0,...,xn1]T的离散余弦变换(DCT)是n维向量y=[y0,...,yn1]T,其中

    y=Cx
    C=2n[121212cosπ2ncos3π2ncos(2n1)π2ncos2π2ncos6π2ncos2(2n1)π2ncos(n1)π2ncos(n1)3π2ncos(n1)(2n1)π2n]

    注意到C是实数正交矩阵,即转置与逆矩阵相同。

     

    定理(DCT插值定理):令x=[x0,...,xn1]T为n维实数向量,定义y=[y0,...,yn1]T=Cx,其中C是n阶DCT变换矩阵,则实数函数

    Pn(t)=1ny0+2nk=1n1ykcosk(2t+1)π2n

    满足Pn(j)=xj,j=0,...,n-1。

     

    举例:用DCT变换计算以下点的插值,(0,1),(1,0),(2,-1),(3,0)。

    先计算4×4的C矩阵,x=[1,0,1,0]T,计算Cx得到y即为插值函数的系数,计算得到的y为

    y=[0.00000.92391.00000.3827]

    所以n=4时得到的插值函数为

    P4(t)=12[0.9239cos(2t+1)π8+cos2(2t+1)π80.3827cos3(2t+1)π8]

    DCT变换和最小二乘近似

    定理(DCT最小二乘近似定理):设x=[x0,...,xn1]T为n维实数向量,令y=[y0,...,yn1]T=Cx,其中C是DCT变换矩阵,则对于任意正整数m≤n,通过选择系数y0,...,ym1,得到的函数

    Pm(t)=1ny0+2nk=1m1ykcosk(2t+1)π2n

    可以使平方误差j=0n1(Pm(j)xj)2最小化。

     

    二维DCT和图像压缩

    二维DCT

    二维DCT变换实际上只是把一维DCT先后应用到二维数据的各个维度上,二维DCT可以用来进行二维网格数据的插值和拟合,计算过程与一维的情形类似。把纵坐标作为第一维,横坐标作为第二维,二维DCT的目的是构造插值函数F(s,t)拟合n2个点(si,tj,xij),其中i,j=0,...,n-1。从最小二乘的观点来看,二维DCT用最优的方法实现了以上的拟合。

    二维DCT是将一维DCT变换先后应用到水平和竖直方向上。考虑二维数据xij组成的矩阵X,先在水平s方向上应用一维DCT变换,首先需要把X转置,再乘以变换矩阵C。结果中的每一列就是X中每行的一维DCT结果,CXT的每一列对应固定的ti,对t方向进行一维DCT变换需要按每行来进行;因此再次转置并乘以C,得到

    C(CXT)T=CXCT

    定义:n×n矩阵X的二维DCT变换(2D-DCT)定义为Y=CXCT

    定义:n×n矩阵Y的二维DCT逆变换是矩阵X=CTYC

    正交变换(例如2D-DCT) 与插值有密切的关系,插值的过程是为了恢复原始数据点,先通过DCT变换得到插值系数,然后用这些系数构造插值函数,再在函数上采样得到数据点。由于C是正交阵,C1=CT,2D-DCT的逆变换可以写成插值的形式,X=CTYC,因为在这个方程中的xij实际上是余弦函数的乘积。

    对于X=CTYC,2D-DCT插值如下

    Pn(i,j)=xij=k=0n1t=0n1CikTyklClj=k=0n1t=0n1CkiyklClj=2nk=0n1t=0n1yklakalcosk(2i+1)π2ncosl(2j+1)π2n

    Pn(i,j)=xij,其中i,j=0,...,n-1。

     

    图像压缩

    一般对8×8的图像块进行DCT变换,将图像中的能量集中到图像块的左上角。

     

    量化

    模q量化

    量化:z=round(yq)

    反量化:y¯=qz

    线性量化定义为矩阵

    qkl=8p(k+l+1)

    其中0≤k,l≤7,其中常量p称为损失参数。

     

    RGB -> GRAY

    Xgray=0.2126R+0.7152G+0.0722B

    RGB -> YUV

    Y=0.299R+0.587G+0.114BU=BYV=RY

    霍夫曼编码

    不过多赘述

     

    改进的DCT和音频压缩

    人的听觉系统对频率非常敏感,在压缩和解压过程中产生的任何瑕疵都能被听出来,基于这个原因,音频压缩中通常会用到一些复杂的技巧来掩盖压缩带来的影响。

    改进的DCT

    定义:离散余弦变换版本4(DCT4)是指对n维向量x=(x0,...,xn1)T,其变换为

    y=Ex

    其中E为n×n矩阵

    Eij=2ncos(i+12)(j+12)πn

    引理:设cj为(扩展的)DCT4矩阵的第j列,则有

    (1)cj=c1j(各列关于j=-1/2对称),

    (2)cj=c2n1j(各列关于j=n-1/2反对称)。

    我们将用DCT4对应的矩阵E来构造MDCT变换,设n为偶数,我们要用列cn2,...,c52n1构造一个新的矩阵。上面的引理说明,对任意整数j,cj都可以用DCT4中的某一列来表示,即0≤i≤n-1中的某个ci,如下图所示。

    image-20231105185249948

    定义:令n为正偶数,向量x=(x0,,x2n1)T修正的DCT(MDCT)为n维向量

    y=Mx

    其中M是n×2n矩阵

    Mij=2ncos(i+12)(j+n2+12)πn

    其中,0≤i≤n-1,0≤j≤2n-1。

    与前面讲到的DCT变换的最大区别是: MDCT把2n维的向量转换为了n维的向量。基于这个原因, MDCT不是可逆的,但是通过重叠这些长度为2n的向量,也可以达到可逆的目的。

    我们可以把MDCT的变换矩阵M用DCT4的列来表示,然后用上面的引理来简化,得到

    M=[cn2c52n1]=[cn2cn1|cnc3n21|c3n2c2n1|c2nc52n1]=[cn2cn1|cn1cn2|cn21c0|c0cn21]

    对n=4的MDCT矩阵

    M=[c2c3|c4c5|c6c7|c8c9]=[c2c3|c3c2|c1c0|c0c1]

    令A和B表示DCT4矩阵的左半部分和右半部分,则E=[A|B]。定义一个置换矩阵R,可以把矩阵的各列顺序颠倒:

    R=[1...1]

    当一个矩阵右乘以矩阵R,则它的各列的顺序会颠倒。若矩阵左乘以R,则其各行的顺序会颠倒。

    利用置换矩阵,可以进一步简化MDCT矩阵的表示

    M=(B|BR|AR|A)

    其中AR和BR是列序颠倒后的A和B矩阵。

    MDCT变换可以用DCT4来表示,令

    x=[x1x2x3x4]

    为一个2n维的向量,其中xi是一个n/2维的向量(注意n为偶数)

    Mx=Bx1Bx2ARx3Ax4=[A|B][Rx3x4x1Rx2]=E[Rx3x4x1Rx2]

    其中E是一个n×n的DCT4变换矩阵,Rx2Rx3表示把x2x3的元素顺序颠倒。这样我们就可以用正交矩阵E来表示M。

    由于MDCT变换的矩阵M是一个n×2n矩阵,而不是方阵,因此它是不可逆的。但是,两个相邻的MDCT 的总阶数为2n,因此若把它们组合起来,就可以完美地重建出输入数据x。

    MDCT的“逆”表示为一个2n×n矩阵N=MT,是M的转置。

    Nij=2ncos(j+12)(i+n2+12)πn
    N=[BTRBTRATAT]

    因为E是正交阵,可以得到

    ATA=IBTB=IATB=BTA=0R2=I
    NM[x1x2x3x4]=[BTRBTRATAT][Bx1Bx2ARx3Ax4]=[x1Rx2Rx1+x2x3+Rx4Rx3+x4]

    在音频压缩中,MDCT变换的向量对应于重叠的数据段。压缩误差会使每段向量之间的连接处产生不连续,而由于向量的长度是固定的,因此这种误差会带来固定频率的噪声。而听觉系统对周期性错误的敏感度要高于视觉系统;实际上,固定周期的错误就是那个频率上的一个声调,而人耳正是感知的各种声调的。如果数据表示采用重叠的方式,令n为正偶数,

    Z1=[x1x2x3x4],Z2=[x3x4x5x6]

    为两个2n维的向量,其中每个xi是一个长度为n/2的向量,向量Z1Z2有一半长度是重叠的,

    NMZ1=[x1Rx2Rx1+x2x3+Rx4Rx3+x4],NMZ2=[x3Rx4Rx3+x4x5+Rx6Rx5+x6]

    NMZ1的下半部分和NMZ2的上半部分求和,可以重建出n维向量[x3x4]:

    [x3x4]=12(NMZ1)n2n1+12(NMZ2)0n1

    位量化

    通过对信号的MDCT进行量化,可以实现音频信号的有损压缩。通过推广图像压缩中使用的量化方法,可以更好地控制信号表示所需的位数。

    首先确定一个实数开区间(-L,L),假设我们要用b位来表示(-L,L)上的一个数,而且可以接受一定程度的误差。我们可以用一个位来表示符号,并用b—1位来对数值进行量化。公式为:

    (-L,L)的b位量化

    量化:z=round(yq),其中q=2L2b1

    反量化:y¯=qz

    特征值与奇异值

    幂迭代方法

    在求解方程中遇到的威尔金森多项式,难以寻找特征值。本节介绍的方法基于对矢量乘上矩阵的高阶幂,这个向量随着幂的升高会变成特征向量,在后面我们将对这个方法进行精化,但是这就是最复杂方法的主要思想。

     

    幂迭代

    幂迭代的动机是与矩阵相乘可以将向量推向主特征向量的方向。

    定义:令A是一个m×m矩阵,A的占优特征值是量级比矩阵A所有其他特征值都大的特征值λ,如果这样的特征值存在,与λ相关特征向量被称为占优特征向量

    幂迭代本质上是一个在每步上进行归一化的不动点迭代。和FPI一样,该方法线性收敛,意味着在收敛过程中,在每个迭代步骤中误差以常数因子降低。在本节的随后部分,将讨论幂迭代具有二阶收敛的变种,该方法称为瑞利商迭代。

    考虑特征值方程xλ=Ax,其中x是特征向量的近似,λ未知,在这种情况下,系数矩阵是n×1矩阵x,法线方程指出最小二乘解是xTxλ=xTAx的解,或者

    λ=xTAxxTx

    这就是瑞利商,给定特征向量近似,瑞利商是特征值的最优近似。

    幂迭代

    给定初始向量x0

    for j = 1,2,3,...

    uj1=xj1/||xj1||2

    xj=Auj1

    λj=uj1TAuj1

    end

    uj=xj/||xj||2

     

    幂迭代的收敛

    定理:令A是一个m×m矩阵,特征值为λ1,...,λm,并满足|λ1|>|λ2||λ3||λm|。假设矩阵A的特征向量张成Rm空间,对于几乎所有的初始向量,幂迭代线性收敛到和λ1相关的特征向量,收敛常数为S=|λ2/λ1|

     

    幂迭代的逆

    幂迭代局限于求解最大(绝对值最大)的特征值,如果幂迭代用于矩阵的逆矩阵,可以找到最小的特征值。

    逆向幂迭代

    给定初始向量x0以及平移s

    for j = 1,2,3,...

    uj1=xj1/||xj1||2

    求解(AsI)xj=uj1

    λj=uj1Txj

    end

    uj=xj/||xj||2

    为了找出矩阵A在实数s附近的特征值,对(AsI)1使用幂迭代得到(AsI)1的最大特征值b。幂迭代可以通过对(AsI)yk+1=xk进行高斯消去得到。则λ=b1+s为矩阵A在s附近的特征值.和λ相关的特征向量可由计算直接给出。

     

    瑞利商迭代

    瑞利商可以同逆向幂迭代同时使用,我们知道它会收敛到和s最接近的特征值对应的特征向量,如果这个距离很小则收敛速度很快。如果在这个过程中任何步骤里,近似的特征值已知,可以使用该特征值作为s,以加速收敛。使用瑞利商作为逆向幂迭代中更新的平移会得到瑞利商迭代(RQI)方法。

    瑞利商迭代

    给定初值x0

    for j=1,2,3,...

    uj1=xj1/||xj1||

    λj1=uj1TAuj1

    求解(Aλj1I)xj=uj1

    end

    uj=xj/||xj||2

     

    QR算法

    本节的目标是推导可以一次找出所有特征值的方法,我们从一个用于对称矩阵的方法开始,然后对该方法进行补充,并用于一般矩阵的特征值求解。对称矩阵容易处理是因为它的特征值为实数,其特征向量构成Rm空间(见附录A)的一组单位正交基。这激发了同时对m个向量使用幂迭代,其中我们需要保持每个向量都与其他向量正交。

    同时迭代

    假设开始有m个两两正交的初始向量,对每个向量使用幂方法一步后得到的向量不再保证彼此正交,所以在每步后重新对m个向量正交化,正交步骤可以看做对结果进行QR分解,如果使用初等基向量作为初始向量,则重新正交化后幂迭代的第一步是AI=Q¯1R1或者

    [A[100]|A[010]||A[001]|]=[q11||qm1][r111r121r1m1r221rmm1]

    qi¯1是幂迭代过程中新的正交单位向量集,其中i=1,...,m。

    归一化同时迭代(NSI)

    Q0¯=I

    for j=1,2,3,...

    AQ¯j=Q¯j+1Rj+1

    end

    在第j步,Q¯j的列是A的特征向量的近似,对角线元素r11j,...,rmmk是近似的特征值。

     

    无移动的QR算法

    RjQj替代Aj

     

    定理:假设A是对称m×m矩阵,特征值为λi,满足|λ1|>|λ2|>>>|λm|。无移动的QR算法可以线性收敛到A的特征值和特征向量。当j→∞时,Aj收敛到对角线矩阵,主对角线上包含所有的特征值,Q¯j=Q1Qj收敛到正交矩阵,对应的列是特征向量。

     

    实数舒尔形式和QR算法

    定义:如果矩阵T为上三角矩阵,且在主对角线上可以出现2×2的块,该矩阵具有实数舒尔形式

    定理:令A是实数元素的方阵,则存在正交矩阵Q以及实数舒尔形式的矩阵T,满足A=QTTQ

    完整的QR算法迭代通过一系列相似变换,将任意矩阵A移动到它的舒尔分解。我们将使用两个步骤,首先将移动的逆向幂迭代思想植入,并加入收缩技术推出移动QR算法,然后推出改进版本允许处理复数特征值。平移版本可以非常直观地写出。每步中都需要使用平移,完成QR分解,然后再平移回来,使用符号表示:

    A0sI=Q1R1A1=R1Q1+sI

    注意到

    A1sI=R1Q1=Q1T(A0sI)Q1=Q1TA0Q1sI

    意味着A1A0相似,因而具有相同的特征值,重复这个步骤,生成一系列Ak矩阵,都与A=A0相似。什么是最优的平移?这带来了特征值计算中的收缩概念。我们将选择平移量为矩阵Ak的右下角元素,这将导致迭代随着矩阵收敛为实数舒尔形式,而将矩阵的最下面一行除了最右边元素之外全部变为0。在这个元素收敛为特征值后,我们去掉矩阵的最后一行和最后一列进行收缩,然后继续寻找其他的特征值。

    原书中有平移QR算法和改进版本的matlab代码

     

    上海森伯格形式

    定义:如果aij=0,i>j+1,m×n矩阵A是上海森伯格形式

    定理:令A是一个方阵,存在正交矩阵Q,满足A=QBQT,并且B是上海森伯格形式。

    可以使用用于构造QR分解的豪斯霍尔德反射构造B,但是在这里有一个主要的差别:现在我们关注在矩阵的左侧和右侧乘上反射子H,这是由于最后想得到一个具有相同的特征值的相似矩阵。正因为如此,我们必须逐步地将0放入矩阵A中。

    原书中给出了对应的matlab代码

     

    奇异值分解

    Rm中对单位球施加m×m矩阵变换得到一个椭球,这个有趣的事实是奇异值分解的基础。对于每个m×n的矩阵A,具有单位正交集{u1,...,um}和{v1,...,vn},以及非负数字s1sn0,满足

    Av1=s1u1Av2=s2u2Avn=snun

    vi被称为矩阵A的右奇异向量,ui是矩阵A的左奇异向量,si是矩阵A的奇异值。

    举例:找出矩阵A的奇异值和奇异向量

    A=[0123000]
    Av1=A[10]=3[010]=s1u1Av2=A[01]=12[100]=s2u2

    在对m×n矩阵A分解的过程中,有一个标准的方式来记录所有的信息,生成一个m×m矩阵U,它的列为左奇异向量ui,一个n×n矩阵V的列为右奇异向量vi,以及一个m×n的对角线矩阵S,其对角线元素是奇异值si,m×n矩阵A的奇异值分解(SVD)

    A=USVT

    找出一般的SVD

    引理:令A是一个m×n矩阵,ATA的特征值非负。

    对于一个m×n的矩阵A,n×n矩阵ATA对称,因而它的特征向量正交,特征值是实数。又特征值是非负实数,因而可以表达为s12≥...≥sn2,其对应的正交特征向量集为{v1,...,vn},使用下面的方向找出ui,其中1≤i≤m:

    如果si0,使用方程siui=Avi定义ui

    如果si=0,选择一个任意的单位向量ui,与u1,...,ui1正交。

    注意ui是两两正交的单位向量(当si=0时,根据定义可得,当si0,有

    (siui)Tsjuj=uiTsiTsjuj=siTsjuiTuj(Avi)TAvj=viTATAvj=λviTvj=0siTsjuiTuj=λviTvj=0

    )因此是Rm的一个正交基。事实上,u1,...,um构成了AAT的正交特征向量。

    定理:令A是m×n矩阵,则存在两个正交基:Rm的{v1,...,vn},与Rm的{u1,...,um},以及实数s1≥...≥sn≥0,满足Avi=siui,1≤i≤min{m,n},V=[v1||vn]的列是右奇异向量,并构成ATA的单位正交向量;U=[u1||um]的列是左奇异向量,构成了ATA的单位正交向量。

    举例:找出2×2矩阵A的奇异值和奇异向量

    A=[0101]
    ATA=[0002]

    的奇异值以降序排列为v1=[0,1]s12=2v2=[1,0]s22=0。奇异值为2与0,根据前面的方向,u1定义为

    2u1=Av1=[11]

    求解u1=[1/2,1/2]T,因为s2为0,所以选择u2u1正交即可,u2=[1/2,1/2]T,SVD为

    [0101]=[2/22/22/22/2][2000][0110]

    注意SVD分解并不唯一,还有另外一种分解

    [0101]=[2/22/22/22/2][2000][0110]

    对于对称矩阵,只需计算自身的特征值和特征向量即可。

     

    SVD的应用

    SVD的性质

    假设A=USVT是奇异值分解,m×n矩阵A的秩是线性无关行的个数(也可以是线性无关列的个数)

    性质1:矩阵A=USVT的秩是S中非0元素的个数。

    性质2:如果A为n×n矩阵,则|det(A)|=s1sn

    性质3:如果A为可逆的m×m矩阵,则A1=VS1UT

    性质4:m×n矩阵A可以写成秩为1的矩阵的和

    A=i=1rsiuiviT

    其中r为A的秩,uivi分别是U和V的第i列。该性质是SVD的低秩近似性质,对于A的最优最小二乘近似为保留的前p项,p为矩阵A的秩,p≤r。

     

    降维

    降维思想是将数据投影到低维空间。假设a1,...,an包含一组m维的向量,在富含数据的应用中,m远小于n,降维的目标是使用n个向量替换a1,...,an,其中新向量的维数p<m,同时最小化该过程中引入的误差。通常我们从一组均值为0的向量开始,如果向量的均值不为0,我们可以减去均值实现,并在后面的阶段再把均值加回去。

    SVD则给出了一种直接完成降维的方式,考虑数据向量为m×n矩阵A=[a1||an],计算奇异值分解A=USVT,令ej表示第j个初等基向量(除了第j个元素为1,其他元素都为0),则Aej=aj,使用性质4的秩近似

    AAp=i=1psiuiviT

    我们可以将aj投影到p维空间,该空间由U的列向量u1,...,up张成

    aj=AejApej

    由于矩阵与ej相乘仅捡出第j列,我们可以这样描述:空间<u1,,up>由左奇异向量u1,...,up张成,这是对于a1,...,an的p维子空间在最小二乘意义上的近似,A的列ai在该空间上的正交投影对应Ap的列。换句话说,一组向量a1,...,an到其最优的最小二乘p维子空间的投影就是矩阵最优的秩p近似矩阵Ap

     

    压缩

    可以利用性质4来压缩矩阵信息,注意到在性质4中的秩1展开中的每一项使用两个向量uivi;以及另一个数字si定义。如果A是一个n×n矩阵,我们可以尝试矩阵A的有损压缩,其中扔掉性质4求和后面的几项,它们具有较小的si。每一项在展开中需要2n+1个数字保存或者传输。

     

    计算SVD

    如果A是一个实数对称矩阵,SVD可以消减为本章前面讨论的特征值计算问题。在这种情况下单位特征向量构成正交基。如果我们定义矩阵V将单位特征向量保存在列上,则AV=US表示特征向量方程,其中S是对角线矩阵,保存特征值的绝对值。对于U的操作和V一样,但是如果特征值是负,则要变换列的符号,由于U和V是正交矩阵,

    A=USVT

    是矩阵A的奇异值分解。

    对于一个一般的非对称m×n矩阵A,确定SVD有两种不同的方法。最明显的方法是构造ATA,但当A的条件数很大时,该方法可能不适合使用。

    考虑如下矩阵

    B=[0ATA0]

    注意,B是对称的(m+n)×(m+n)矩阵,因而,它具有实数特征值和作为基的特征向量。令[v,w]表示一个(m+n)向量,对应B的特征向量,则

    [ATwAv]=[0ATA0][vw]=λ[vw]

    或者Av=λw,左侧乘上AT得到

    ATAv=λATw=λ2v

    表明w是ATA的特征向量,对应的特征值为λ2,注意到我们可以确定ATA的特征值和特征向量,但是不需要构造ATA

    因而,另一种更好的计算奇异值和奇异向量的方法是,将对称矩阵B转化为上海森伯格形式。由于对称,上海森伯格形式的矩阵和三对角线矩阵等价,则该方法与平移QR算法一样可以找出特征值,这是奇异值的平方,特征向量中n个最大的项是奇异向量v、尽管这种方式看起来使得矩阵的大小加倍,但是它避免对条件数的不必要提高。

     

    最优化

    最优化指的是找出实数函数的极大或者极小值,该函数称为目标函数。由于定位函数f(z)的极大值与找出函数-f(x)的极小值等价,在推导计算方式时仅考虑最小化问题就足够了。

     

    不使用导数的无约束优化

     

    黄金分割搜索

    定义:当区间[a,b]上只有一个极大或者极小值,并且f在其他点上严格升高或降低,连续函数f(x)被称为区间[a,b]上的单峰函数

    定理:从初始区间[a,b]开始,黄金分割搜索k步之后,最后区间中点和最小值之间的差异为gk(ba)/2,其中g=(51)/20.618

    黄金分割搜索

    给出单峰函数f,在区间[a,b]上具有极小值

    for i = 1,2,3,...

    g=(51)/2

    if f(a+(1-g)(b-a)) < f(a+g(b-a))

    b = a+g(b-a)

    else

    a = a+(1-g)(b-a)

    end

    end

    最终区间[a,b]中包含极小值

     

    持续的抛物线插值

    持续抛物线插值

    从近似极小值r,s,t开始 for i = 1,2,3,... x=r+s2(f(s)f(r))(tr)(ts)2[(sr)(f(t)f(s))(f(s)f(r))(ts)] t = s s = r r = x end

    x是对极小值的新的近似。在SPI中,新的x可以替换r,s,t中离当前最远,或者最差的一个点,并重复进行该步骤。对于SPI不能保证收敛,这与黄金分割搜索不同。但是,如果收敛的话,由于该方法更好地使用了函数求值的信息,所以速度会更快。

     

    Nelder-Mead搜索

    Nelder-Mead搜索试图将一个多面体滚到一个尽可能低的水平。由于这个原因,它也被称为单纯型下山法。它没有使用目标函数的导数信息。假设需要最小化的函数f具有n个变量,该方法首先需要n+1个属于Rm的初始估计向量x1,...,xn+1,这些点构成n维的单纯型。测试单纯型的顶点,并根据它们的值以升序排列y1<y2<...<yn+1=yh。最差的单纯型向量xh=xn+1根据下面的流程图进行替换。首先我们定义略去xh的单纯型平面的重心x¯。然后我们测试在反射点xr=2x¯xh上的函数值yr=f(xr),如图(a)所示,如果新的值yry1<y2<yn范围中,我们使用xr替换最差的点xn,使用函数值对顶点排序,并重复前面的步骤。

    image-20231106101117408

    具体代码见原书P503页

     

    使用导数的无约束优化

     

    牛顿方法

    如果函数连续可微,可以对导数求导,则优化问题可以表示为求解根的问题。

    在一个连续可微函数f(x)的局部极小值点x处,一阶导数必然为0。可以使用求解方程组的方法求解f'(x)=0。如果目标函数单峰,在区间中具有极小值,则使用最小值x附近的初始估计开始牛顿方法的计算,这将会收敛到x。对f'(x)=0应用牛顿方法得到如下迭代

    xk+1=xkf(xk)f(xk)

    当遇到一个多变量函数时,可以使用多变量的牛顿方法,设置导数为0,有

    f=0

    其中

    f=[fx1(x1,,xn),,fxn(x1,,xn)]

    表示梯度。

    使用方程组一章中的牛顿方法可以求解这个问题,设F(x)=f(x),牛顿方法的迭代步中将设置xk+1=xk+v,其中v是DF(xk)v=F(xk)的解。梯度的雅可比矩阵DF为

    Hf=DF=[2fx1x22fx1xn2fxnx12fxnxn]

    这是函数f的海森矩阵,因而牛顿步骤为

    {Hf(xk)v=f(xk)xk+1=xk+v

    最速下降

    最速下降,也称为梯度搜索,背后的基本想法是将当前点在在最速下降的方向上移动以找出函数最小值。由于f指向f的最速生长方向,反方向f就是最速下降方向。

    最速下降

    for i = 0,1,2,...

    v=f(xi)

    对于标量s=s最小化f(xsv)

    xi+1=xisv

    end

     

    共轭梯度搜索

    当A为对称正定矩阵时,求解Ax=w等价于找出抛物面的极小值,如在二维中,线性方程组

    [abbc][x1x2]=[ef]

    的解是如下抛物面的极小值

    f(x1,x2)=12ax12+bx1x2+12cx22ex1fx2

    这是由于f的梯度可以表示为线性方程组

    [abbc][x1x2][ef]=[ax1+bx2ebx1+cx2f]=f=[fx1fx2]

    在极小值处的梯度为0,正定意味着抛物面的凹面向上。

    一个关键的观测时线性方程组的余项r=w-Ax为f(x),这是函数f在x点的最速下降方向。假设以及选择了一个搜索方向,使用向量d标记。为了在那个方向找到最小值,需要找出α使得函数h(α)=f(x+αd)最小。将导数设置为0以找出最小:

    0=fd=(A(x+αd)(e,f)T)d=(αAdr)Td

    这意味着

    α=rTddTAd=rTrdTAd

    共轭梯度搜索

    x0为初始估计,设d0=r0=f

    for i = 1,2,3,...

    αi=α使得f(xi1+αdi1)最小

    xi=xi1+αidi1

    ri=f(xi)

    βi=riTriri1Tri1

    di=ri+βidi1

    end

     

    附录

    分块乘法

    分块计算A×B,分块的唯一要求是:对A矩阵的各列的分组方法要与B矩阵各行的分组方法一致,如

    image-20231106114503654